画像処理装置、画像処理装置の信頼性評価方法および画像処理プログラム
【課題】推定オペレータを用いた分光特性の推定における観測ノイズの影響を予測すること。
【解決手段】影響度算出部１４３は、推定オペレータからノイズに対する相対的な影響度を算出する。また、重み算出部１４５は、ノイズに対する相対的な影響度をもとに色素量推定に用いる重み係数を算出する。そして、スペクトル推定部１４７は、対象標本画像の画素値に基づいて対象標本のスペクトル（分光透過率）を推定し、色素量推定部１４９は、重み係数を用い、対象標本の染色に用いたヘマトキシリンおよびエオジンの各色素の基準分光特性をもとに対象標本の色素量を推定する。
【解決手段】影響度算出部１４３は、推定オペレータからノイズに対する相対的な影響度を算出する。また、重み算出部１４５は、ノイズに対する相対的な影響度をもとに色素量推定に用いる重み係数を算出する。そして、スペクトル推定部１４７は、対象標本画像の画素値に基づいて対象標本のスペクトル（分光透過率）を推定し、色素量推定部１４９は、重み係数を用い、対象標本の染色に用いたヘマトキシリンおよびエオジンの各色素の基準分光特性をもとに対象標本の色素量を推定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、推定オペレータを用いて分光特性の推定を行う画像処理装置、この画像処理装置の信頼性評価方法および画像処理プログラムに関する。
【背景技術】
【０００２】
被写体に固有の物理的性質を表す物理量の一つに分光透過率スペクトルがある。分光透過率は、各波長における入射光に対する透過光の割合を表す物理量であり、ＲＧＢ値等の照明光の変化に依存する色情報とは異なり、外因的影響によって値が変化しない物体固有の情報である。このため、分光透過率は、被写体自体の色を再現するための情報として様々な分野で利用されている。例えば、生体組織標本、特に病理標本を用いた病理診断の分野では、標本を撮像した画像の解析に分光特性の一例として分光透過率の推定技術が利用されている。
【０００３】
病理診断では、臓器摘出によって得たブロック標本や針生検によって得た病理標本を厚さ数ミクロン程度に薄切した後、様々な所見を得るために顕微鏡を用いて拡大観察することが広く行われている。中でも光学顕微鏡を用いた透過観察は、機材が比較的安価で取り扱いが容易である上、歴史的に古くから行われてきたこともあって、最も普及している観察方法の一つである。この場合、薄切された標本は光を殆ど吸収及び散乱せず無色透明に近いため、観察に先立って色素による染色を施すのが一般的である。
【０００４】
染色手法としては種々のものが提案されており、その総数は１００種類以上にも達するが、特に病理標本に関しては、色素として青紫色のヘマトキシリンと赤色のエオジンの２つを用いるヘマトキシリン−エオジン染色（以下、「Ｈ＆Ｅ染色」と称す。）が標準的に用いられている。
【０００５】
ヘマトキシリンは植物から採取された天然の物質であり、それ自身には染色性はない。しかし、その酸化物であるヘマチンは好塩基性の色素であり、負に帯電した物質と結合する。細胞核に含まれるデオキシリボ核酸（ＤＮＡ）は、構成要素として含むリン酸基によって負に帯電しているため、ヘマチンと結合して青紫色に染色される。なお、前述の通り、染色性を有するのはヘマトキシリンでは無く、その酸化物であるヘマチンであるが、色素の名称としてはヘマトキシリンを用いるのが一般的であるため、以下それに従う。一方エオジンは、好酸性の色素であり、正に帯電した物質と結合する。アミノ酸やタンパク質が正負どちらに帯電するかはｐＨ環境に影響を受け、酸性下では正に帯電する傾向が強くなる。このため、エオジン溶液に酢酸を加えて用いることがある。細胞質に含まれるタンパク質は、エオジンと結合して赤から薄赤に染色される。
【０００６】
Ｈ＆Ｅ染色後の標本（染色標本）では、細胞核や骨組織等が青紫色に、細胞質や結合組織、赤血球等が赤色に染色され、容易に視認できるようになる。この結果、観察者は、細胞核等の組織を構成する要素の大きさや位置関係等を把握でき、染色標本の状態を形態学的に判断することが可能となる。
【０００７】
染色標本の観察は、観察者の目視によるものの他、この染色標本をマルチバンド撮像して外部装置の表示画面に表示することによっても行われている。表示画面に表示する場合には、撮像したマルチバンド画像から標本各点の分光透過率を推定する処理や、推定した分光透過率に基づいて標本を染色している色素の色素量を推定する処理、推定した色素量に基づいて画像の色を補正する処理等が行われ、カメラの特性や染色状態のばらつき等が補正されて、表示用の標本のＲＧＢ画像が合成される。図１２は、合成されたＲＧＢ画像の一例を示す図である。色素量の推定を適切に行えば、濃く染色された標本や薄く染色された標本を、適切に染色された標本と同等の色を有する画像に補正することができる。したがって、染色標本の分光透過率を高精度に推定することが、染色標本に固定された色素量の推定や、染色ばらつきの補正等の高精度化に繋がる。
【０００８】
標本のマルチバンド画像から標本各点の分光透過率を推定する手法としては、例えば、主成分分析による推定法（例えば、非特許文献１参照）や、ウィナー（Wiener）推定による推定法（例えば、非特許文献２参照）等が挙げられる。ウィナー推定は、ノイズの重畳された観測信号から原信号を推定する線形フィルタ手法の一つとして広く知られており、観測対象の統計的性質とノイズ（観測ノイズ）の特性とを考慮して誤差の最小化を行う手法である。カメラからの信号には何らかのノイズが含まれるため、ウィナー推定は原信号を推定する手法として極めて有用である。
【０００９】
ここで、ウィナー推定法によって標本のマルチバンド画像から標本各点の分光透過率を推定する方法について説明する。
【００１０】
先ず、標本のマルチバンド画像を撮像する。例えば、特許文献１に開示されている技術を用い、１６枚のバンドパスフィルタをフィルタホイールで回転させて切り替えながら、面順次方式でマルチバンド画像を撮像する。これにより、標本の各点において１６バンドの画素値を有するマルチバンド画像が得られる。なお、色素は、本来観察対象となる染色標本内に３次元的に分布しているが、通常の透過観察系ではそのまま３次元像として捉えることはできず、標本内を透過した照明光をカメラの撮像素子上に投影した２次元像として観察される。したがって、ここでいう各点は、投影された撮像素子の各画素に対応する標本上の点を意味している。
【００１１】
撮像されたマルチバンド画像の任意の点ｘについて、バンドｂにおける画素値ｇ（ｘ，ｂ）と、対応する標本上の点の分光透過率ｔ（ｘ，λ）との間には、カメラの応答システムに基づく次式（１）の関係が成り立つ。
【数１】
λは波長、ｆ（ｂ，λ）はｂ番目のフィルタの分光透過率、ｓ（λ）はカメラの分光感度特性、ｅ（λ）は照明の分光放射特性、ｎ（ｂ）はバンドｂにおける観測ノイズをそれぞれ表す。ｂはバンドを識別する通し番号であり、ここでは１≦ｂ≦１６を満たす整数値である。
【００１２】
実際の計算では、式（１）を波長方向に離散化した次式（２）を用いる。
Ｇ（ｘ）＝ＦＳＥＴ（ｘ）＋Ｎ ・・・（２）
波長方向のサンプル点数をＤ、バンド数をＢとすれば（ここではＢ＝１６）、Ｇ（ｘ）は、点ｘにおける画素値ｇ（ｘ，ｂ）に対応するＢ行１列の行列である。同様に、Ｔ（ｘ）は、ｔ（ｘ，λ）に対応するＤ行１列の行列、Ｆは、ｆ（ｂ，λ）に対応するＢ行Ｄ列の行列である。一方、Ｓは、Ｄ行Ｄ列の対角行列であり、対角要素がｓ（λ）に対応している。同様に、Ｅは、Ｄ行Ｄ列の対角行列であり、対角要素がｅ（λ）に対応している。Ｎは、ｎ（ｂ）に対応するＢ行１列の行列である。なお、式（２）では、行列を用いて複数のバンドに関する式を集約しているため、バンドを表す変数ｂが陽に記述されていない。また、波長λに関する積分は行列の積に置き換えられている。
【００１３】
ここで、表記を簡単にするため、次式（３）で定義される行列Ｈを導入する。Ｈはシステム行列とも呼ばれる。
Ｈ＝ＦＳＥ ・・・（３）
【００１４】
次に、ウィナー推定を用いて、撮像したマルチバンド画像から標本各点における分光透過率を推定する。分光透過率の推定値Ｔ＾（ｘ）は、次式（４）で計算することができる。なお、Ｔ＾は、Ｔの上に推定値を表すハット（＾）が付いていることを示す。
【数２】
【００１５】
ここで、Ｗは次式（５）で表され、「ウィナー推定行列」あるいは「ウィナー推定に用いる推定オペレータ」と呼ばれる。以下の説明では、Ｗを単に「推定オペレータ」と称す。
【数３】
ＲSSは、Ｄ行Ｄ列の行列であり、標本の分光透過率の自己相関行列を表す。また、ＲNNは、Ｂ行Ｂ列の行列であり、撮像に使用するカメラのノイズの自己相関行列を表す。
【００１６】
このようにして分光透過率Ｔ＾（ｘ）を推定したならば、次に、このＴ＾（ｘ）に基づいて対応する標本点における色素量を推定する。推定の対象とする色素は、第１の色素量に相当する細胞核を染色したヘマトキシリン、第２の色素量に相当する細胞質を染色したエオジン、第３の色素量に相当する赤血球を染色したエオジンの３種類である。ここで、ヘマトキシリンを色素Ｈ、細胞質を染色したエオジンを色素Ｅ、赤血球を染色したエオジンを色素Ｒと略記する。なお、厳密には、染色を施さない状態であっても赤血球はそれ自身特有の色を有しており、Ｈ＆Ｅ染色後は、赤血球自身の色と染色過程において変化したエオジンの色が重畳して観察される。このため、正確には両者を併せたものを色素Ｒと呼称する。
【００１７】
一般に、光を透過する物質では、波長λ毎の入射光の強度Ｉ0（λ）と射出光の強度Ｉ(λ)との間に、次式（６）で表されるLambert-Beer則が成り立つことが知られている。
【数４】
ｋ（λ）は波長に依存して決まる物質固有の値、ｄは物質の厚さをそれぞれ表す。また、式（６）の左辺は分光透過率を意味している。
【００１８】
Ｈ＆Ｅ染色された対象標本が、色素Ｈ、色素Ｅ、色素Ｒの３種類の色素で染色されている場合、Lambert-Beer則により各波長λにおいて次式（７）が成立する。
【数５】
ここでｋH（λ），ｋE（λ），ｋR（λ）は、それぞれ色素Ｈ、色素Ｅ、色素Ｒに対応したｋ（λ）を表す。またｄH，ｄE，ｄRは、マルチバンド画像の各画像位置に対応する標本各点における色素Ｈ、色素Ｅ、色素Ｒの仮想的な厚さを表す。本来色素は、標本中に分散して存在するため、厚さという概念は正確ではないが、標本が単一の色素で染色されていると仮定した場合と比較して、どの程度の量の色素が存在しているかを表す相対的な色素量の指標となる。すなわち、ｄH，ｄE，ｄRはそれぞれ色素Ｈ、色素Ｅ、色素Ｒの色素量を表しているといえる。なお、ｋH（λ），ｋE（λ），ｋR（λ）は、単一の色素で染色した標本を予め用意し、その分光透過率を分光計で測定することによって、Lambert-Beer則から容易に求めることができる。
【００１９】
式（７）の両辺の対数を取ると、次式（８）となる。
【数６】
【００２０】
式（４）を用いて推定された分光透過率データＴ＾（ｘ）の波長λに対応する要素をｔ＾（ｘ，λ）とし、これを式（８）に代入すると、次式（９）を得る。
【数７】
ここで、式（９）において未知変数はｄH，ｄE，ｄRの３つであるから、少なくとも３つの異なる波長λに対して式（９）を連立させれば、これらを解くことができる。より精度を高めるために、４つ以上の異なる波長λに対して式（９）を連立させ、重回帰分析を行ってもよい。
【００２１】
そして、このようにして標本点における色素量を推定したならば、推定した色素量をもとに、標本各点における色素量を適正な染色状態に調整し、標本の画像を合成する。
【００２２】
【特許文献１】特開平７−１２０３２４号公報
【非特許文献１】“Development of support systems for pathology using spectral transmittance - The quantification method of stain conditions”，Proceedings of SPIE，Vol.4684，2002，p.1516-1523
【非特許文献２】“Color Correction of Pathological Images Based on Dye Amount Quantification”，OPTICAL REVIEW，Vol.12，No.4，2005，p.293-300
【発明の開示】
【発明が解決しようとする課題】
【００２３】
非特許文献１に開示されている手法に従い、推定オペレータを用いてマルチバンド画像から標本各点の分光透過率を推定する場合、式（１）に示したカメラの応答システムに基づく関係が成り立ち、点ｘのバンドｂにおける画素値ｇ（ｘ，ｂ）は、フィルタの分光透過率ｆ（ｂ，λ）、カメラの分光感度特性ｓ（λ）、照明の分光放射特性ｅ（λ）および分光透過率ｔ（ｘ，λ）によって求まる信号値と、観測ノイズｎ（ｂ）とによって定まる。しかしながら、従来の分光特性（分光透過率）の推定におけるウィナー推定は、ノイズの影響を考慮した推定方法であるが、算出したウィナー推定行列がノイズにどのように影響を受けるかは十分には解明されておらず、分光特性の推定精度が低下する恐れがあり、色素量推定の精度を確保できない場合があった。また、この結果、推定オペレータを用いた分光特性の推定を行う装置の信頼性を適正に判断できないという問題が生じていた。また、推定に必要なバンド数を求めるには実験を繰り返す必要があった。
【００２４】
本発明は、上記した従来の問題点に鑑みて為されたものであり、推定オペレータを用いた分光特性の推定におけるノイズの影響度を予測することができる画像処理装置、この画像処理装置の信頼性評価方法および画像処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【００２５】
上記した課題を解決し、目的を達成するため、本発明にかかる画像処理装置は、推定オペレータを用いて分光特性の推定を行う画像処理装置において、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析手段を備えることを特徴とする。
【００２６】
また、本発明にかかる画像処理装置は、染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定手段と、分光特性の波長毎の重み係数を取得する重み係数取得手段と、前記分光特性推定手段によって推定された分光特性と前記重み係数取得手段によって取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定手段と、を備えることを特徴とする。
【００２７】
また、本発明にかかる画像処理装置の信頼性評価方法は、推定オペレータを用いて分光特性の推定を行う画像処理装置の信頼性評価方法であって、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析し、解析したノイズに対する相対的な影響度をもとに前記画像処理装置の装置信頼性を評価することを特徴とする。
【００２８】
また、本発明にかかる画像処理プログラムは、推定オペレータを用いて分光特性の推定を行うコンピュータに、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析ステップを実行させることを特徴とする。
【００２９】
また、本発明にかかる画像処理プログラムは、コンピュータに、染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定ステップと、分光特性の波長毎の重み係数を取得する重み係数取得ステップと、前記分光特性推定ステップで推定された分光特性と前記重み係数取得ステップで取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定ステップと、を実行させることを特徴とする。
【発明の効果】
【００３０】
本発明によれば、分光特性の推定に用いる推定オペレータからノイズに対する相対的な影響度を解析することができるので、推定オペレータを用いた分光特性の推定におけるノイズの影響度を予測することができるという効果を奏する。
【発明を実施するための最良の形態】
【００３１】
以下、図面を参照し、本発明の好適な実施の形態について詳細に説明する。本実施の形態では、Ｈ＆Ｅ染色された生体組織標本（染色標本）を撮像対象とし、撮像したマルチバンド画像をもとに、ウィナー推定を用いて染色標本の標本各点のスペクトル推定を行い、標本各点の色素量の推定を行う。なお、この実施の形態によって本発明が限定されるものではない。また、図面の記載において、同一部分には同一の符号を付して示している。
【００３２】
（実施の形態１）
図１は、実施の形態１の画像処理装置の構成を説明する模式図である。図１に示すように、画像処理装置１は、パソコン等のコンピュータで構成され、染色標本のマルチバンド画像を取得する画像取得部１１０を備える。
【００３３】
画像取得部１１０は、画像取得動作を行ってＨ＆Ｅ染色された色素量の推定対象の染色標本（以下、「対象標本」という。）を撮像し、６バンドのマルチバンド画像を取得する。この画像取得部１１０は、ＣＣＤ等の撮像素子等を備えたＲＧＢカメラ１１１、対象標本Ｓが載置される標本保持部１１３、標本保持部１１３上の対象標本Ｓを透過照明する照明部１１５、対象標本Ｓからの透過光を集光して結像させる光学系１１７、結像する光の波長帯域を所定範囲に制限するためのフィルタ部１１９等を備える。
【００３４】
ＲＧＢカメラ１１１は、デジタルカメラ等で広く用いられているものであり、モノクロの撮像素子上にモザイク状にＲＧＢのカラーフィルタを配置したものである。このＲＧＢカメラ１１１は、撮像される画像の中心が照明光の光軸上に位置するように設置される。図２は、カラーフィルタの配列例およびＲＧＢ各バンドの画素配列を模式的に示す図である。この場合、各画素はＲ，Ｇ，Ｂいずれかの成分しか撮像することはできないが、近傍の画素値を利用することで、不足するＲ，Ｇ，Ｂ成分が補間される。この手法は、例えば特許第３５１００３７号公報で開示されている。なお、３ＣＣＤタイプのカメラを使用すれば、最初から各画素におけるＲ，Ｇ，Ｂ成分を取得できる。実施の形態１では、いずれの撮像方式を用いても構わないが、以下ではＲＧＢカメラ１１１で撮像された画像の各画素においてＲ，Ｇ，Ｂ成分が取得できているものとする。
【００３５】
フィルタ部１１９は、それぞれ異なる分光透過率特性を有する２枚の光学フィルタ１１９１ａ，１１９１ｂを具備しており、これらが回転式の光学フィルタ切替部１１９３に保持されて構成されている。図３−１は、一方の光学フィルタ１１９１ａの分光透過率特性を示す図であり、図３−２は、他方の光学フィルタ１１９１ｂの分光透過率特性を示す図である。例えば先ず、光学フィルタ１１９１ａを用いて第１の撮像を行う。次いで、光学フィルタ切替部１１９３の回転によって使用する光学フィルタを光学フィルタ１１９１ｂに切り替え、光学フィルタ１１９１ｂを用いて第２の撮像を行う。この第１の撮像及び第２の撮像によって、それぞれ３バンドの画像が得られ、両者の結果を合わることによって６バンドのマルチバンド画像が得られる。なお、光学フィルタの数は２枚に限定されるものではなく、３枚以上の光学フィルタを用いることができる。取得された染色標本のマルチバンド画像である染色標本画像は、対象標本画像として画像処理装置１の記憶部１５０に保持される。
【００３６】
この画像取得部１１０において、照明部１１５によって照射された照明光は、標本保持部１１３上に載置された対象標本Ｓを透過する。そして、対象標本Ｓを透過した透過光は、光学系１１７及び光学フィルタ１１９１ａ，１１９１ｂを経由した後、ＲＧＢカメラ１１１の撮像素子上に結像する。光学フィルタ１１９１ａ，１１９１ｂを具備するフィルタ部１１９は、照明部１１５からＲＧＢカメラ１１１に至る光路上のいずれかの位置に設置されていればよい。照明部１１５からの照明光を、光学系１１７を介してＲＧＢカメラ１１１で撮像する際の、Ｒ，Ｇ，Ｂ各バンドの分光感度の例を、図４に示す。
【００３７】
図５は、実施の形態１の画像処理装置１の機能構成を説明するブロック図である。実施の形態１では、画像処理装置１は、図１に示して説明した画像取得部１１０と、入力部１２０と、表示部１３０と、演算部１４０と、記憶部１５０と、装置各部を制御する制御部１６０とを備える。
【００３８】
制御部１６０は、ＣＰＵ等のハードウェアによって実現される。この制御部１６０は、入力部１２０から入力される操作信号や画像取得部１１０から入力される画像データ、記憶部１５０に格納されるプログラムやデータ等に基づいて画像処理装置１を構成する各部への指示やデータの転送等を行い、画像処理装置１全体の動作を統括的に制御する。また、制御部１６０は、画像取得部１１０の動作を制御して対象標本画像を取得するマルチバンド画像取得制御部１６１を含む。
【００３９】
入力部１２０は、例えば、キーボードやマウス、タッチパネル、各種スイッチ等の各種入力装置によって実現されるものであり、操作入力に応じた操作信号を制御部１６０に出力する。表示部１３０は、ＬＣＤやＥＬＤ等の表示装置によって実現されるものであり、制御部１６０から入力される表示信号に基づいて各種画面を表示する。
【００４０】
演算部１４０は、ＣＰＵ等のハードウェアによって実現される。この演算部１４０は、推定オペレータを算出する推定オペレータ算出部１４１と、推定オペレータからノイズに対する相対的な影響度を算出する影響度解析手段としての影響度算出部１４３と、ノイズに対する相対的な影響度をもとに色素量推定に用いる重み係数を算出する重み係数取得手段および重み係数設定手段としての重み算出部１４５と、対象標本画像の画素値に基づいて対象標本のスペクトル（分光透過率）を推定する分光特性推定手段としてのスペクトル推定部１４７と、対象標本の染色に用いたヘマトキシリンおよびエオジンの各色素の基準分光特性をもとに対象標本の色素量を推定する色素量推定手段としての色素量推定部１４９とを含む。
【００４１】
記憶部１５０は、更新記憶可能なフラッシュメモリ等のＲＯＭやＲＡＭといった各種ＩＣメモリ、内蔵或いはデータ通信端子で接続されたハードディスク、ＣＤ−ＲＯＭ等の情報記憶媒体及びその読取装置等によって実現されるものであり、画像処理装置１の動作にかかるプログラムや、画像処理装置１の備える種々の機能を実現するためのプログラム、これらプログラムの実行にかかるデータ等が格納される。例えば、対象標本画像の画像データや、推定オペレータ算出部１４１によって算出された推定オペレータ（ウィナー推定行列）のデータ、重み算出部１４５によって算出された重み係数のデータ等が格納される。また、推定オペレータからノイズに対する相対的な影響度を解析して色素量推定に用いる重み係数を算出する処理（以下、「重み係数算出処理」と称す。）や、重み係数算出処理で算出された重み係数を用いて対象標本の色素量を推定する処理（以下、「重み付き色素量推定処理」と称す。）を実現するための画像処理プログラム１５１が格納される。
【００４２】
次に、実施の形態１の画像処理装置１が行う処理手順について説明する。実施の形態１の画像処理装置１は、重み係数算出処理および重み付き色素量推定処理を実行する。なお、ここで説明する処理は、記憶部１５０に格納された画像処理プログラム１５１に従って画像処理装置１の各部が動作することによって実現される。
【００４３】
先ず、重み係数算出処理について、図６を参照して説明する。図６に示すように、重み係数算出処理では、先ず、推定オペレータ算出部１４１が、推定オペレータＷを算出する（ステップＳ１１）。具体的には、背景技術で示した次式（５）に従って推定オペレータＷを算出する。
【数８】
【００４４】
ここで、背景技術で示したように、次式（３）で定義されるシステム行列Ｈを導入する。
Ｈ＝ＦＳＥ ・・・（３）
光学フィルタ１１９１ａ，１１９１ｂの分光透過率Ｆ、ＲＧＢカメラ１１１の分光感度特性Ｓおよび単位時間当たりの照明の分光放射特性Ｅ（＾）は、使用する機器を選定の後、分光計等を用いて予め測定しておく。なお、ここでは、光学系１１７の分光透過率は１.０と近似しているが、この近似値１.０からの乖離が許容できない場合には、光学系１１７の分光透過率も予め測定し、照明の分光放射特性Ｅに乗じればよい。また、標本の分光透過率の自己相関行列ＲSSおよびＲＧＢカメラ１１１のノイズの自己相関行列ＲNNについても、事前に測定しておく。ＲSSは、Ｈ＆Ｅ染色された典型的な標本を用意し、分光計によって複数の点の分光透過率を測定して自己相関行列を求めることによって得られる。ＲNNは、標本無しの状態で画像取得部１１０によってマルチバンド画像を取得し、得られた６バンドのマルチバンド画像の各バンドについて画素値の分散を求め、これを対角成分とする行列を生成することによって得られる。ただし、バンド間でノイズの相関はないと仮定している。
【００４５】
図７は、横軸を波長λ（ｎｍ）とし、縦軸を波長λにおけるバンドｂの推定オペレータ値ｗλbとして、推定オペレータＷの各波長λに対応する要素の値（推定オペレータ値ｗλb）をバンド毎（バンド１〜バンド６）にプロットして推定オペレータＷをグラフ化した図である。算出された推定オペレータＷのデータは、記憶部１５０に保持される。
【００４６】
次に、図６に示すように、影響度算出部１４３が、ステップＳ１１で算出された推定オペレータＷからバンド毎のノイズに対する相対的な影響度を算出する（ステップＳ１３）。具体的には、次式（１０）に従ってバンド毎のノイズに対する相対的な影響度を算出する。
【数９】
Ｂはバンド数を表す。ｃλbは、推定オペレータ値ｗλbのノイズに対する相対的な影響度である。ｃλbの値は、波長λにおけるバンドｂの推定オペレータ値ｗλbを、波長λにおける全バンドの推定オペレータ値の合計値で割ることによって得た値（波長λにおける全バンドの推定オペレータ値に対する所定のバンドｂにおける推定オペレータ値ｗλbの相対的な大きさ）であり、すなわち、推定オペレータＷを用いて推定される染色標本の分光透過率の波長λにおける推定値に対する推定オペレータ値ｗλbの相対的な影響度を意味する。ここで、推定オペレータ値ｗλbは、入力値が原信号でも観測ノイズでも同様に影響する。したがって、ｃλbの値は、バンド毎の観測ノイズに対する相対的な影響度を意味する。
【００４７】
図８は、横軸を波長λ（ｎｍ）、縦軸をバンド毎のノイズに対する相対的な影響度ｃλbとして、算出されたバンド毎のノイズに対する相対的な影響度ｃλbの一例をグラフ化した図である。図８に示す例では、５４０〜５６０ｎｍの波長帯域におけるバンド毎のノイズに対する相対的な影響度ｃλbの値が、他の波長帯域と比較して高い値を示している。このように、バンド毎のノイズに対する相対的な影響度ｃλbは、ノイズ（観測ノイズ）に対する波長依存性を示しており、このバンド毎のノイズに対する相対的な影響度ｃλbを算出して解析することによって、推定オペレータＷを用いたスペクトル推定における観測ノイズの影響度、すなわち推定オペレータＷのノイズ感度を予測できる。具体的には、その値の大きい波長帯域では観測ノイズによる影響を受け易く、観測ノイズに対するノイズ耐性が低いと予測できる。一方、バンド毎のノイズに対する相対的な影響度ｃλbが小さい波長帯域では観測ノイズによる影響を受け難く、観測ノイズに対するノイズ耐性が高いと予測できる。すなわち、図８の例では、５４０〜５６０ｎｍの波長帯域が、他の波長帯域と比較して観測ノイズによる影響を受け易いと予測できる。実施の形態１では、このバンド毎のノイズに対する相対的な影響度ｃλbを用いて観測ノイズによる影響を受け難い波長帯域におけるスペクトル推定の結果が色素量推定に反映され易く、観測ノイズによる影響を受け易い波長帯域におけるスペクトル推定の結果が色素量推定に反映され難くなるように重み係数ωλを波長毎に算出する。
【００４８】
すなわち、図６に示すように、重み算出部１４５が、対象標本の色素量を推定する際に用いる重み係数を波長毎に算出する（ステップＳ１５）。具体的には先ず、次式（１１）に従って、バンド毎のノイズに対する相対的な影響度ｃλbから波長毎のノイズに対する相対的な影響度ｃλを算出する。
【数１０】
なお、ここでは、先鋭的特徴を捉えるため、バンド毎のノイズに対する相対的な影響度ｃλbの絶対値の最大値を用いたが、絶対値の平均値やノルムを用いて波長毎のノイズに対する相対的な影響度ｃλを求めることとしてもよい。例えば、図８に示す例においてバンド２に着目すれば、５２０〜５４０ｎｍの波長帯域で得られたプロットＰ１１に対応するノイズに対する相対的な影響度ｃλbがノイズに対する相対的な影響度ｃλとして算出される。また、バンド５に着目すれば、５４０〜５６０ｎｍの波長帯域で得られたプロットＰ１３に対応するノイズに対する相対的な影響度ｃλbがノイズに対する相対的な影響度ｃλとして算出される。
【００４９】
次いで、波長λにおける重み係数ωλを算出する。ここで、色素量推定の際に、ノイズに対する相対的な影響度ｃλの低い波長帯域に設定される重み係数ωλが大きく、ノイズに対する相対的な影響度の高い波長帯域に設定される重み係数ωλが小さくなるように波長λにおける重み係数ωλを算出する。例えば、次式（１２），（１３），（１４）に従って波長λにおける重み係数ωλを算出し、波長λにおけるノイズに対する相対的な影響度ｃλが低い場合には重み係数ωλの値を大きくし、ノイズに対する相対的な影響度ｃλが高い場合には重み係数ωλの値を小さくする。
【数１１】
【００５０】
また、αは所定の閾値であり、このαの値によって、ノイズに対する相対的な影響度ｃλの高い波長帯域を色素量の推定にどこまで利用するかを決定する。αの値を高くしすぎると重み係数ωλが常に１となり、重み係数ωλの役割を果たさなくなる。一方、αの値を低くしすぎるとｃλ＞１の波長帯域全てに重み係数ωλが設定されてしまう。例えば、図８に示した例では、ほとんどの波長帯域に重み係数ωλが設定されることとなる。したがって、予めノイズに対する相対的な影響度を考慮せずに色素量推定を行った結果から決定する等して最も推定精度が高まるようにαの値を設定する必要がある。
【００５１】
なお、ここでは、ノイズに対する相対的な影響度ｃλが相対的に高い波長帯域の重み係数ωλに０以上１以下の値を設定することとしたが、０を設定することによって、ノイズに対する相対的な影響度ｃλが相対的に高い波長帯域を色素量の推定に用いないこととしてもよい。また、ノイズに対する相対的な影響度ｃλが相対的に高い波長帯域にのみ重み係数ωλを設定することとしたが、全ての波長帯域に対して相対的に重み係数ωλを設定してもよい。算出された重み係数ω1，ω2，・・・，ωDは、記憶部１５０に保持される。
【００５２】
そして、次式（１５）に従って、推定オペレータＷのノイズに対する相対的な影響度ｃを算出する。
【数１２】
なお、ここでは、先鋭的特徴を捉えるため、ｃ1，ｃ2，・・・，ｃDの最大値を用いたが、平均値やノルムを用いて推定オペレータＷのノイズに対する相対的な影響度ｃを求めることとしてもよい。例えば、図８に示す例では、５４０〜５６０ｎｍの波長帯域で得られたプロットＰ１３に対応するバンド５のノイズに対する相対的な影響度ｃλが、ノイズに対する相対的な影響度ｃとして算出される。
【００５３】
ここで求めたノイズに対する相対的な影響度ｃは、推定オペレータＷの信頼性を表す尺度として用いることができ、例えば画像処理装置１の装置信頼性を評価するために参照される。すなわち、ノイズに対する相対的な影響度ｃの値が小さい程スペクトル推定における観測ノイズの影響度は低く、推定オペレータＷの信頼性が高い。一方、ノイズに対する相対的な影響度ｃの値が大きい場合には、スペクトル推定における観測ノイズの影響度が高くなり、推定オペレータＷの信頼性は低下する。このノイズに対する相対的な影響度ｃを用いて画像処理装置１の装置信頼性を評価する際には、その値が小さければ画像処理装置１の装置信頼性を高いと評価する。一方、ノイズに対する相対的な影響度ｃの値が大きければ画像処理装置１の装置信頼性は低いと評価する。
【００５４】
続いて、重み付き色素量推定処理の手順について、図９を参照して説明する。図９に示すように、重み付き色素量推定処理では、先ず、マルチバンド画像取得制御部１６１が、画像取得部１１０の動作を制御して色素量の推定対象の対象標本をマルチバンド撮像し、対象標本画像を取得する（ステップＳ２１）。
【００５５】
次に、スペクトル推定部１４７が、ステップＳ２１で取得した対象標本画像の画素値をもとに、対象標本のスペクトル（分光透過率）を推定する（ステップＳ２３）。具体的には、図６のステップＳ１１で算出した推定オペレータＷを用いる。そして、背景技術で示した次式（４）に従い、対象標本画像の推定対象画素である任意の点ｘにおける画素の画素値の行列表現Ｇ（ｘ）から、対応する対象標本の標本点における分光透過率の推定値Ｔ＾（ｘ）を推定する。得られた分光透過率の推定値Ｔ＾（ｘ）は、記憶部１５０に格納される。
【数１３】
【００５６】
そして、色素量推定部１４９が、ステップＳ２３で推定した分光透過率の推定値Ｔ＾（ｘ）に基づき、重み係数算出処理（図６）で算出した重み係数ωλを用いて対象標本の色素量を推定する（ステップＳ２５）。ここで、色素量推定部１４９は、対象標本の染色に用いたヘマトキシリンおよびエオジンの各色素の基準分光特性をもとに、対象標本画像の任意の点ｘに対応する標本点における細胞核を染色したヘマトキシリン（色素Ｈ）の色素量、細胞質を染色したエオジン（色素Ｅ）の色素量、および赤血球を染色したエオジン（色素Ｒ）の色素量を推定する。具体的には、対象標本画像の点ｘにおける分光透過率の推定値Ｔ＾（ｘ）に基づいて、点ｘに対応する対象標本の標本点に固定された色素Ｈ，色素Ｅ，色素Ｒそれぞれの色素量を推定する。すなわち、背景技術で示した次式（９）を複数の波長λに関して連立させ、ｄH，ｄE，ｄRについて解く。
【数１４】
【００５７】
例として、３つの波長λ1，λ2，λ3について式（９）を連立させた場合を考えると、次式（１６）のように行列表記できる。
【数１５】
ここで、式（１６）を次式（１７）に置き換える。
【数１６】
この関係式（１７）から、最小二乗法を用いて色素量ｄH，ｄE，ｄRを算出する。最小二乗法とは、単回帰式において誤差の二乗和を最小にするようにｄ（ｘ）を決定する方法であり、一般的には次式（１８）で算出できる。
【数１７】
【００５８】
実施の形態１では、重み係数算出処理で算出した重み係数ωλ（ω1，ω2，・・・，ωD）を加味して色素量ｄH，ｄE，ｄRを算出する。このため、次式（１９），（２０）を用いる。
【数１８】
ωは、重み係数ωλに対応するＤ行Ｄ列の行列であり、ｄｉａｇ（）は、対角行列を表す。以上のようにして重み係数ωλを用いて推定された対象標本画像の点ｘに対応する染色標本の標本点における色素量ｄH，ｄE，ｄRは、記憶部１５０に保持される。そして、例えば、ここで推定された色素量ｄH，ｄE，ｄRに基づいて対象標本画像の色が補正され、カメラの特性や染色状態のばらつき等が補正されて表示用のＲＧＢ画像が合成される。このＲＧＢ画像は、表示部１３０に画面表示されて病理診断に利用される。
【００５９】
実施の形態１によれば、推定オペレータＷからバンド毎のノイズに対する相対的な影響度ｃλbを算出し、波長毎のノイズに対する相対的な影響度ｃλを算出することができ、推定オペレータＷを用いたスペクトル推定における観測ノイズの影響度、すなわち推定オペレータＷのノイズ感度を予測できる。また、波長毎のノイズに対する相対的な影響度ｃλから波長毎に重み係数ωλを算出し、算出した重み係数ωλを用いて対象標本の色素量を推定することができる。具体的には、ノイズに対する相対的な影響度ｃλの低い波長帯域の重み係数ωλを大きく、ノイズに対する相対的な影響度ｃλの高い波長帯域の重み係数ωλを小さくして色素量の推定を行うことができるので、推定オペレータＷを用いたスペクトル推定時に受ける観測ノイズの影響を加味して色素量を推定することができ、推定精度の低下を軽減できる。また、ノイズに対する相対的な影響度ｃλをもとに算出したノイズに対する相対的な影響度ｃを用いて画像処理装置１の装置信頼性を評価することができ、推定オペレータＷを用いたスペクトル推定時に受ける観測ノイズの影響下における装置の信頼性を適正に判断することができる。
【００６０】
なお、実施の形態１では、画像処理装置１において、推定オペレータ算出部１４１が推定オペレータＷを算出し、重み算出部１４５が重み係数ωλを算出することによって重み係数ωλを取得する場合について説明したが、推定オペレータＷや重み係数ωλは、推定の都度毎回算出しなくてもよい。すなわち、予め推定オペレータＷおよび／または重み係数ωλを算出し、記憶部１５０に格納しておく構成としてもよい。そして、スペクトル推定の際には記憶部１５０から推定オペレータＷのデータを読み出して用い、色素量推定の際には記憶部１５０から重み係数ωλのデータを読み出して用いる構成としてもよい。
【００６１】
（実施の形態２）
次に、実施の形態２について説明する。図１０は、実施の形態２の画像処理装置１ｂの機能構成を説明するブロック図である。なお、実施の形態１で説明した構成と同一の構成については、同一の符号を付する。実施の形態２では、画像処理装置１ｂは、図１に示して説明した画像取得部１１０と、入力部１２０と、表示部１３０と、演算部１４０ｂと、記憶部１５０ｂと、装置各部を制御する制御部１６０とを備える。そして、演算部１４０ｂは、影響度算出部１４３ｂと、重み算出部１４５ｂと、スペクトル推定部１４７ｂと、色素量推定部１４９ｂとを含む。
【００６２】
また、記憶部１５０ｂには、予め求められた推定オペレータＷ（ウィナー推定行列）のデータが推定オペレータデータ１５３ｂとして格納されるとともに、推定オペレータデータ１５３ｂに格納された推定オペレータＷと対象標本をマルチバンド撮像した対象標本画像の画素値とからノイズに対する相対的な影響度を解析して色素量推定に用いる重み係数を算出し、この重み係数を用いて対象標本の色素量を推定する処理を実現するための画像処理プログラム１５１ｂが格納される。
【００６３】
図１１は、実施の形態２の画像処理装置１ｂが行う処理手順を示すフローチャートである。なお、ここで説明する処理は、記憶部１５０ｂに格納された画像処理プログラム１５１ｂに従って画像処理装置１ｂの各部が動作することによって実現される。
【００６４】
先ず、マルチバンド画像取得制御部１６１が、画像取得部１１０の動作を制御して色素量の推定対象の対象標本をマルチバンド撮像し、対象標本画像を取得する（ステップＳ３１）。次に、スペクトル推定部１４７ｂが、推定オペレータデータ１５３ｂから推定オペレータＷを読み出し、この推定オペレータＷを用い、ステップＳ３１で取得した対象標本画像の画素値をもとに対象標本のスペクトル（分光透過率）を推定する（ステップＳ３３）。
【００６５】
続いて、影響度算出部１４３ｂが、推定オペレータデータ１５３ｂから読み出した推定オペレータＷと、ステップＳ３１で取得した対象標本画像の画素値とからバンド毎のノイズに対する相対的な影響度を算出する（ステップＳ３５）。具体的には、次式（２１）に従ってバンド毎のノイズに対する相対的な影響度ｃλbを算出する。
【数１９】
ｇ(ｂ)は、バンドｂの対象標本画像の推定対象画素である点ｘの画素値である。このようにして実施の形態２では、推定オペレータ値ｗλbに加えて、推定対象画素の画素値を用いてバンド毎のノイズに対する相対的な影響度ｃλbを算出する。
【００６６】
そして、重み算出部１４５ｂが、波長毎の重み係数を算出する（ステップＳ３７）。重み係数の算出の手法は実施の形態１と同様であり、ステップＳ３５で求めたバンド毎のノイズに対する相対的な影響度ｃλbから式（１１）に従って波長毎のノイズに対する相対的な影響度ｃλを算出し、波長毎のノイズに対する相対的な影響度ｃλから式（１２）〜（１４）に従って波長λにおける重み係数ωλを算出する。また、式（１５）に従って、画像処理装置１ｂの装置信頼性を評価するための値である推定オペレータＷのノイズに対する相対的な影響度ｃを算出する。ここで求めたノイズに対する相対的な影響度ｃは、実施の形態１と同様にして画像処理装置１ｂの装置信頼性を評価するために参照される。
【００６７】
そして、色素量推定部１４９ｂが、ステップＳ３３で推定した分光透過率の推定値Ｔ＾（ｘ）に基づき、ステップＳ３７で算出した重み係数ωλを用いて対象標本の色素量を推定する（ステップＳ３９）。色素量の推定手法は、実施の形態１と同様である。
【００６８】
実施の形態２によれば、実施の形態１と同様の効果を奏するとともに、推定オペレータＷとバンド毎の対象標本画像の画素値ｇ(ｂ)とからバンド毎のノイズに対する相対的な影響度ｃλbを算出し、波長毎のノイズに対する相対的な影響度ｃλを算出することができ、推定オペレータＷを用いたスペクトル推定における観測ノイズの影響、すなわち推定オペレータＷのノイズ感度を対象標本画像の画素値を用いて予測することができる。
【００６９】
なお、実施の形態２では、推定対象画素の画素値を用いてバンド毎のノイズに対する相対的な影響度ｃλbを算出する場合について説明したが、例えば、予め複数の画素値についてバンド毎のノイズに対する相対的な影響度ｃλbを解析し、重み係数ωλを算出して記憶部１５０ｂに格納しておくこととしてもよい。そして、色素量の推定時には、推定対象画素の画素値に応じて最適な重み係数ωλを記憶部１５０ｂから読み出して用いることとしてもよい。例えば、所定幅で画素値をサンプリングし、サンプリングした各画素値についてそれぞれ重み係数ωλを算出しておくこととしてもよい。そして、色素量の推定時には、推定対象画素の画素値と最も値が近い画素値の重み係数ωλを用いて色素量の推定を行う。あるいは、対象標本画像に映る細胞核や細胞質、赤血球、背景といった組織毎に重み係数ωλを算出しておくこととしてもよい。具体的には、組織に応じた代表的な画素値を用いて組織毎の重み係数ωλとして算出しておくこととしてもよい。そして、色素量の推定時には、推定対象画素の画素値から対応する標本点の組織を特定し、特定した組織に応じた重み係数ωλを用いて色素量の推定を行う。
【００７０】
また、実施の形態２では、推定パラーメータＷを予め算出しておく場合について説明したが、推定の都度算出することとしてもよい。
【００７１】
また、上記の実施の形態では、病理標本を撮像したマルチバンド画像から分光透過率のスペクトル特徴値を推定する場合について説明したが、分光特性値として、分光反射率のスペクトル特徴値を推定する場合にも、同様に適用できる。
【図面の簡単な説明】
【００７２】
【図１】画像処理装置の構成を示す図である。
【図２】カラーフィルタの配列例およびＲＧＢ各バンドの画素配列を模式的に示す図である。
【図３−１】一の光学フィルタの分光透過率特性を示す図である。
【図３−２】他の光学フィルタの分光透過率特性を示す図である。
【図４】Ｒ，Ｇ，Ｂ各バンドの分光感度の例を示す図である。
【図５】実施の形態１の画像処理装置の機能構成を示すブロック図である。
【図６】重み係数算出処理の処理手順を示すフローチャートである。
【図７】推定オペレータをグラフ化した図である。
【図８】バンド毎のノイズに対する相対的な影響度の一例をグラフ化した図である。
【図９】重み付き色素量推定処理の処理手順を示すフローチャートである。
【図１０】実施の形態２の画像処理装置の機能構成を示すブロック図である。
【図１１】実施の形態２の画像処理装置が行う処理手順を示すフローチャートである。
【図１２】ＲＧＢ画像の一例を示す図である。
【符号の説明】
【００７３】
１ 画像処理装置
１１０ 画像取得部
１２０ 入力部
１３０ 表示部
１４０ 演算部
１４１ 推定オペレータ算出部
１４３ 影響度算出部
１４５ 重み算出部
１４７ スペクトル推定部
１４９ 色素量推定部
１５０ 記憶部
１５１ 画像処理プログラム
１６０ 制御部
１６１ マルチバンド画像取得制御部
【技術分野】
【０００１】
本発明は、推定オペレータを用いて分光特性の推定を行う画像処理装置、この画像処理装置の信頼性評価方法および画像処理プログラムに関する。
【背景技術】
【０００２】
被写体に固有の物理的性質を表す物理量の一つに分光透過率スペクトルがある。分光透過率は、各波長における入射光に対する透過光の割合を表す物理量であり、ＲＧＢ値等の照明光の変化に依存する色情報とは異なり、外因的影響によって値が変化しない物体固有の情報である。このため、分光透過率は、被写体自体の色を再現するための情報として様々な分野で利用されている。例えば、生体組織標本、特に病理標本を用いた病理診断の分野では、標本を撮像した画像の解析に分光特性の一例として分光透過率の推定技術が利用されている。
【０００３】
病理診断では、臓器摘出によって得たブロック標本や針生検によって得た病理標本を厚さ数ミクロン程度に薄切した後、様々な所見を得るために顕微鏡を用いて拡大観察することが広く行われている。中でも光学顕微鏡を用いた透過観察は、機材が比較的安価で取り扱いが容易である上、歴史的に古くから行われてきたこともあって、最も普及している観察方法の一つである。この場合、薄切された標本は光を殆ど吸収及び散乱せず無色透明に近いため、観察に先立って色素による染色を施すのが一般的である。
【０００４】
染色手法としては種々のものが提案されており、その総数は１００種類以上にも達するが、特に病理標本に関しては、色素として青紫色のヘマトキシリンと赤色のエオジンの２つを用いるヘマトキシリン−エオジン染色（以下、「Ｈ＆Ｅ染色」と称す。）が標準的に用いられている。
【０００５】
ヘマトキシリンは植物から採取された天然の物質であり、それ自身には染色性はない。しかし、その酸化物であるヘマチンは好塩基性の色素であり、負に帯電した物質と結合する。細胞核に含まれるデオキシリボ核酸（ＤＮＡ）は、構成要素として含むリン酸基によって負に帯電しているため、ヘマチンと結合して青紫色に染色される。なお、前述の通り、染色性を有するのはヘマトキシリンでは無く、その酸化物であるヘマチンであるが、色素の名称としてはヘマトキシリンを用いるのが一般的であるため、以下それに従う。一方エオジンは、好酸性の色素であり、正に帯電した物質と結合する。アミノ酸やタンパク質が正負どちらに帯電するかはｐＨ環境に影響を受け、酸性下では正に帯電する傾向が強くなる。このため、エオジン溶液に酢酸を加えて用いることがある。細胞質に含まれるタンパク質は、エオジンと結合して赤から薄赤に染色される。
【０００６】
Ｈ＆Ｅ染色後の標本（染色標本）では、細胞核や骨組織等が青紫色に、細胞質や結合組織、赤血球等が赤色に染色され、容易に視認できるようになる。この結果、観察者は、細胞核等の組織を構成する要素の大きさや位置関係等を把握でき、染色標本の状態を形態学的に判断することが可能となる。
【０００７】
染色標本の観察は、観察者の目視によるものの他、この染色標本をマルチバンド撮像して外部装置の表示画面に表示することによっても行われている。表示画面に表示する場合には、撮像したマルチバンド画像から標本各点の分光透過率を推定する処理や、推定した分光透過率に基づいて標本を染色している色素の色素量を推定する処理、推定した色素量に基づいて画像の色を補正する処理等が行われ、カメラの特性や染色状態のばらつき等が補正されて、表示用の標本のＲＧＢ画像が合成される。図１２は、合成されたＲＧＢ画像の一例を示す図である。色素量の推定を適切に行えば、濃く染色された標本や薄く染色された標本を、適切に染色された標本と同等の色を有する画像に補正することができる。したがって、染色標本の分光透過率を高精度に推定することが、染色標本に固定された色素量の推定や、染色ばらつきの補正等の高精度化に繋がる。
【０００８】
標本のマルチバンド画像から標本各点の分光透過率を推定する手法としては、例えば、主成分分析による推定法（例えば、非特許文献１参照）や、ウィナー（Wiener）推定による推定法（例えば、非特許文献２参照）等が挙げられる。ウィナー推定は、ノイズの重畳された観測信号から原信号を推定する線形フィルタ手法の一つとして広く知られており、観測対象の統計的性質とノイズ（観測ノイズ）の特性とを考慮して誤差の最小化を行う手法である。カメラからの信号には何らかのノイズが含まれるため、ウィナー推定は原信号を推定する手法として極めて有用である。
【０００９】
ここで、ウィナー推定法によって標本のマルチバンド画像から標本各点の分光透過率を推定する方法について説明する。
【００１０】
先ず、標本のマルチバンド画像を撮像する。例えば、特許文献１に開示されている技術を用い、１６枚のバンドパスフィルタをフィルタホイールで回転させて切り替えながら、面順次方式でマルチバンド画像を撮像する。これにより、標本の各点において１６バンドの画素値を有するマルチバンド画像が得られる。なお、色素は、本来観察対象となる染色標本内に３次元的に分布しているが、通常の透過観察系ではそのまま３次元像として捉えることはできず、標本内を透過した照明光をカメラの撮像素子上に投影した２次元像として観察される。したがって、ここでいう各点は、投影された撮像素子の各画素に対応する標本上の点を意味している。
【００１１】
撮像されたマルチバンド画像の任意の点ｘについて、バンドｂにおける画素値ｇ（ｘ，ｂ）と、対応する標本上の点の分光透過率ｔ（ｘ，λ）との間には、カメラの応答システムに基づく次式（１）の関係が成り立つ。
【数１】
λは波長、ｆ（ｂ，λ）はｂ番目のフィルタの分光透過率、ｓ（λ）はカメラの分光感度特性、ｅ（λ）は照明の分光放射特性、ｎ（ｂ）はバンドｂにおける観測ノイズをそれぞれ表す。ｂはバンドを識別する通し番号であり、ここでは１≦ｂ≦１６を満たす整数値である。
【００１２】
実際の計算では、式（１）を波長方向に離散化した次式（２）を用いる。
Ｇ（ｘ）＝ＦＳＥＴ（ｘ）＋Ｎ ・・・（２）
波長方向のサンプル点数をＤ、バンド数をＢとすれば（ここではＢ＝１６）、Ｇ（ｘ）は、点ｘにおける画素値ｇ（ｘ，ｂ）に対応するＢ行１列の行列である。同様に、Ｔ（ｘ）は、ｔ（ｘ，λ）に対応するＤ行１列の行列、Ｆは、ｆ（ｂ，λ）に対応するＢ行Ｄ列の行列である。一方、Ｓは、Ｄ行Ｄ列の対角行列であり、対角要素がｓ（λ）に対応している。同様に、Ｅは、Ｄ行Ｄ列の対角行列であり、対角要素がｅ（λ）に対応している。Ｎは、ｎ（ｂ）に対応するＢ行１列の行列である。なお、式（２）では、行列を用いて複数のバンドに関する式を集約しているため、バンドを表す変数ｂが陽に記述されていない。また、波長λに関する積分は行列の積に置き換えられている。
【００１３】
ここで、表記を簡単にするため、次式（３）で定義される行列Ｈを導入する。Ｈはシステム行列とも呼ばれる。
Ｈ＝ＦＳＥ ・・・（３）
【００１４】
次に、ウィナー推定を用いて、撮像したマルチバンド画像から標本各点における分光透過率を推定する。分光透過率の推定値Ｔ＾（ｘ）は、次式（４）で計算することができる。なお、Ｔ＾は、Ｔの上に推定値を表すハット（＾）が付いていることを示す。
【数２】
【００１５】
ここで、Ｗは次式（５）で表され、「ウィナー推定行列」あるいは「ウィナー推定に用いる推定オペレータ」と呼ばれる。以下の説明では、Ｗを単に「推定オペレータ」と称す。
【数３】
ＲSSは、Ｄ行Ｄ列の行列であり、標本の分光透過率の自己相関行列を表す。また、ＲNNは、Ｂ行Ｂ列の行列であり、撮像に使用するカメラのノイズの自己相関行列を表す。
【００１６】
このようにして分光透過率Ｔ＾（ｘ）を推定したならば、次に、このＴ＾（ｘ）に基づいて対応する標本点における色素量を推定する。推定の対象とする色素は、第１の色素量に相当する細胞核を染色したヘマトキシリン、第２の色素量に相当する細胞質を染色したエオジン、第３の色素量に相当する赤血球を染色したエオジンの３種類である。ここで、ヘマトキシリンを色素Ｈ、細胞質を染色したエオジンを色素Ｅ、赤血球を染色したエオジンを色素Ｒと略記する。なお、厳密には、染色を施さない状態であっても赤血球はそれ自身特有の色を有しており、Ｈ＆Ｅ染色後は、赤血球自身の色と染色過程において変化したエオジンの色が重畳して観察される。このため、正確には両者を併せたものを色素Ｒと呼称する。
【００１７】
一般に、光を透過する物質では、波長λ毎の入射光の強度Ｉ0（λ）と射出光の強度Ｉ(λ)との間に、次式（６）で表されるLambert-Beer則が成り立つことが知られている。
【数４】
ｋ（λ）は波長に依存して決まる物質固有の値、ｄは物質の厚さをそれぞれ表す。また、式（６）の左辺は分光透過率を意味している。
【００１８】
Ｈ＆Ｅ染色された対象標本が、色素Ｈ、色素Ｅ、色素Ｒの３種類の色素で染色されている場合、Lambert-Beer則により各波長λにおいて次式（７）が成立する。
【数５】
ここでｋH（λ），ｋE（λ），ｋR（λ）は、それぞれ色素Ｈ、色素Ｅ、色素Ｒに対応したｋ（λ）を表す。またｄH，ｄE，ｄRは、マルチバンド画像の各画像位置に対応する標本各点における色素Ｈ、色素Ｅ、色素Ｒの仮想的な厚さを表す。本来色素は、標本中に分散して存在するため、厚さという概念は正確ではないが、標本が単一の色素で染色されていると仮定した場合と比較して、どの程度の量の色素が存在しているかを表す相対的な色素量の指標となる。すなわち、ｄH，ｄE，ｄRはそれぞれ色素Ｈ、色素Ｅ、色素Ｒの色素量を表しているといえる。なお、ｋH（λ），ｋE（λ），ｋR（λ）は、単一の色素で染色した標本を予め用意し、その分光透過率を分光計で測定することによって、Lambert-Beer則から容易に求めることができる。
【００１９】
式（７）の両辺の対数を取ると、次式（８）となる。
【数６】
【００２０】
式（４）を用いて推定された分光透過率データＴ＾（ｘ）の波長λに対応する要素をｔ＾（ｘ，λ）とし、これを式（８）に代入すると、次式（９）を得る。
【数７】
ここで、式（９）において未知変数はｄH，ｄE，ｄRの３つであるから、少なくとも３つの異なる波長λに対して式（９）を連立させれば、これらを解くことができる。より精度を高めるために、４つ以上の異なる波長λに対して式（９）を連立させ、重回帰分析を行ってもよい。
【００２１】
そして、このようにして標本点における色素量を推定したならば、推定した色素量をもとに、標本各点における色素量を適正な染色状態に調整し、標本の画像を合成する。
【００２２】
【特許文献１】特開平７−１２０３２４号公報
【非特許文献１】“Development of support systems for pathology using spectral transmittance - The quantification method of stain conditions”，Proceedings of SPIE，Vol.4684，2002，p.1516-1523
【非特許文献２】“Color Correction of Pathological Images Based on Dye Amount Quantification”，OPTICAL REVIEW，Vol.12，No.4，2005，p.293-300
【発明の開示】
【発明が解決しようとする課題】
【００２３】
非特許文献１に開示されている手法に従い、推定オペレータを用いてマルチバンド画像から標本各点の分光透過率を推定する場合、式（１）に示したカメラの応答システムに基づく関係が成り立ち、点ｘのバンドｂにおける画素値ｇ（ｘ，ｂ）は、フィルタの分光透過率ｆ（ｂ，λ）、カメラの分光感度特性ｓ（λ）、照明の分光放射特性ｅ（λ）および分光透過率ｔ（ｘ，λ）によって求まる信号値と、観測ノイズｎ（ｂ）とによって定まる。しかしながら、従来の分光特性（分光透過率）の推定におけるウィナー推定は、ノイズの影響を考慮した推定方法であるが、算出したウィナー推定行列がノイズにどのように影響を受けるかは十分には解明されておらず、分光特性の推定精度が低下する恐れがあり、色素量推定の精度を確保できない場合があった。また、この結果、推定オペレータを用いた分光特性の推定を行う装置の信頼性を適正に判断できないという問題が生じていた。また、推定に必要なバンド数を求めるには実験を繰り返す必要があった。
【００２４】
本発明は、上記した従来の問題点に鑑みて為されたものであり、推定オペレータを用いた分光特性の推定におけるノイズの影響度を予測することができる画像処理装置、この画像処理装置の信頼性評価方法および画像処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【００２５】
上記した課題を解決し、目的を達成するため、本発明にかかる画像処理装置は、推定オペレータを用いて分光特性の推定を行う画像処理装置において、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析手段を備えることを特徴とする。
【００２６】
また、本発明にかかる画像処理装置は、染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定手段と、分光特性の波長毎の重み係数を取得する重み係数取得手段と、前記分光特性推定手段によって推定された分光特性と前記重み係数取得手段によって取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定手段と、を備えることを特徴とする。
【００２７】
また、本発明にかかる画像処理装置の信頼性評価方法は、推定オペレータを用いて分光特性の推定を行う画像処理装置の信頼性評価方法であって、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析し、解析したノイズに対する相対的な影響度をもとに前記画像処理装置の装置信頼性を評価することを特徴とする。
【００２８】
また、本発明にかかる画像処理プログラムは、推定オペレータを用いて分光特性の推定を行うコンピュータに、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析ステップを実行させることを特徴とする。
【００２９】
また、本発明にかかる画像処理プログラムは、コンピュータに、染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定ステップと、分光特性の波長毎の重み係数を取得する重み係数取得ステップと、前記分光特性推定ステップで推定された分光特性と前記重み係数取得ステップで取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定ステップと、を実行させることを特徴とする。
【発明の効果】
【００３０】
本発明によれば、分光特性の推定に用いる推定オペレータからノイズに対する相対的な影響度を解析することができるので、推定オペレータを用いた分光特性の推定におけるノイズの影響度を予測することができるという効果を奏する。
【発明を実施するための最良の形態】
【００３１】
以下、図面を参照し、本発明の好適な実施の形態について詳細に説明する。本実施の形態では、Ｈ＆Ｅ染色された生体組織標本（染色標本）を撮像対象とし、撮像したマルチバンド画像をもとに、ウィナー推定を用いて染色標本の標本各点のスペクトル推定を行い、標本各点の色素量の推定を行う。なお、この実施の形態によって本発明が限定されるものではない。また、図面の記載において、同一部分には同一の符号を付して示している。
【００３２】
（実施の形態１）
図１は、実施の形態１の画像処理装置の構成を説明する模式図である。図１に示すように、画像処理装置１は、パソコン等のコンピュータで構成され、染色標本のマルチバンド画像を取得する画像取得部１１０を備える。
【００３３】
画像取得部１１０は、画像取得動作を行ってＨ＆Ｅ染色された色素量の推定対象の染色標本（以下、「対象標本」という。）を撮像し、６バンドのマルチバンド画像を取得する。この画像取得部１１０は、ＣＣＤ等の撮像素子等を備えたＲＧＢカメラ１１１、対象標本Ｓが載置される標本保持部１１３、標本保持部１１３上の対象標本Ｓを透過照明する照明部１１５、対象標本Ｓからの透過光を集光して結像させる光学系１１７、結像する光の波長帯域を所定範囲に制限するためのフィルタ部１１９等を備える。
【００３４】
ＲＧＢカメラ１１１は、デジタルカメラ等で広く用いられているものであり、モノクロの撮像素子上にモザイク状にＲＧＢのカラーフィルタを配置したものである。このＲＧＢカメラ１１１は、撮像される画像の中心が照明光の光軸上に位置するように設置される。図２は、カラーフィルタの配列例およびＲＧＢ各バンドの画素配列を模式的に示す図である。この場合、各画素はＲ，Ｇ，Ｂいずれかの成分しか撮像することはできないが、近傍の画素値を利用することで、不足するＲ，Ｇ，Ｂ成分が補間される。この手法は、例えば特許第３５１００３７号公報で開示されている。なお、３ＣＣＤタイプのカメラを使用すれば、最初から各画素におけるＲ，Ｇ，Ｂ成分を取得できる。実施の形態１では、いずれの撮像方式を用いても構わないが、以下ではＲＧＢカメラ１１１で撮像された画像の各画素においてＲ，Ｇ，Ｂ成分が取得できているものとする。
【００３５】
フィルタ部１１９は、それぞれ異なる分光透過率特性を有する２枚の光学フィルタ１１９１ａ，１１９１ｂを具備しており、これらが回転式の光学フィルタ切替部１１９３に保持されて構成されている。図３−１は、一方の光学フィルタ１１９１ａの分光透過率特性を示す図であり、図３−２は、他方の光学フィルタ１１９１ｂの分光透過率特性を示す図である。例えば先ず、光学フィルタ１１９１ａを用いて第１の撮像を行う。次いで、光学フィルタ切替部１１９３の回転によって使用する光学フィルタを光学フィルタ１１９１ｂに切り替え、光学フィルタ１１９１ｂを用いて第２の撮像を行う。この第１の撮像及び第２の撮像によって、それぞれ３バンドの画像が得られ、両者の結果を合わることによって６バンドのマルチバンド画像が得られる。なお、光学フィルタの数は２枚に限定されるものではなく、３枚以上の光学フィルタを用いることができる。取得された染色標本のマルチバンド画像である染色標本画像は、対象標本画像として画像処理装置１の記憶部１５０に保持される。
【００３６】
この画像取得部１１０において、照明部１１５によって照射された照明光は、標本保持部１１３上に載置された対象標本Ｓを透過する。そして、対象標本Ｓを透過した透過光は、光学系１１７及び光学フィルタ１１９１ａ，１１９１ｂを経由した後、ＲＧＢカメラ１１１の撮像素子上に結像する。光学フィルタ１１９１ａ，１１９１ｂを具備するフィルタ部１１９は、照明部１１５からＲＧＢカメラ１１１に至る光路上のいずれかの位置に設置されていればよい。照明部１１５からの照明光を、光学系１１７を介してＲＧＢカメラ１１１で撮像する際の、Ｒ，Ｇ，Ｂ各バンドの分光感度の例を、図４に示す。
【００３７】
図５は、実施の形態１の画像処理装置１の機能構成を説明するブロック図である。実施の形態１では、画像処理装置１は、図１に示して説明した画像取得部１１０と、入力部１２０と、表示部１３０と、演算部１４０と、記憶部１５０と、装置各部を制御する制御部１６０とを備える。
【００３８】
制御部１６０は、ＣＰＵ等のハードウェアによって実現される。この制御部１６０は、入力部１２０から入力される操作信号や画像取得部１１０から入力される画像データ、記憶部１５０に格納されるプログラムやデータ等に基づいて画像処理装置１を構成する各部への指示やデータの転送等を行い、画像処理装置１全体の動作を統括的に制御する。また、制御部１６０は、画像取得部１１０の動作を制御して対象標本画像を取得するマルチバンド画像取得制御部１６１を含む。
【００３９】
入力部１２０は、例えば、キーボードやマウス、タッチパネル、各種スイッチ等の各種入力装置によって実現されるものであり、操作入力に応じた操作信号を制御部１６０に出力する。表示部１３０は、ＬＣＤやＥＬＤ等の表示装置によって実現されるものであり、制御部１６０から入力される表示信号に基づいて各種画面を表示する。
【００４０】
演算部１４０は、ＣＰＵ等のハードウェアによって実現される。この演算部１４０は、推定オペレータを算出する推定オペレータ算出部１４１と、推定オペレータからノイズに対する相対的な影響度を算出する影響度解析手段としての影響度算出部１４３と、ノイズに対する相対的な影響度をもとに色素量推定に用いる重み係数を算出する重み係数取得手段および重み係数設定手段としての重み算出部１４５と、対象標本画像の画素値に基づいて対象標本のスペクトル（分光透過率）を推定する分光特性推定手段としてのスペクトル推定部１４７と、対象標本の染色に用いたヘマトキシリンおよびエオジンの各色素の基準分光特性をもとに対象標本の色素量を推定する色素量推定手段としての色素量推定部１４９とを含む。
【００４１】
記憶部１５０は、更新記憶可能なフラッシュメモリ等のＲＯＭやＲＡＭといった各種ＩＣメモリ、内蔵或いはデータ通信端子で接続されたハードディスク、ＣＤ−ＲＯＭ等の情報記憶媒体及びその読取装置等によって実現されるものであり、画像処理装置１の動作にかかるプログラムや、画像処理装置１の備える種々の機能を実現するためのプログラム、これらプログラムの実行にかかるデータ等が格納される。例えば、対象標本画像の画像データや、推定オペレータ算出部１４１によって算出された推定オペレータ（ウィナー推定行列）のデータ、重み算出部１４５によって算出された重み係数のデータ等が格納される。また、推定オペレータからノイズに対する相対的な影響度を解析して色素量推定に用いる重み係数を算出する処理（以下、「重み係数算出処理」と称す。）や、重み係数算出処理で算出された重み係数を用いて対象標本の色素量を推定する処理（以下、「重み付き色素量推定処理」と称す。）を実現するための画像処理プログラム１５１が格納される。
【００４２】
次に、実施の形態１の画像処理装置１が行う処理手順について説明する。実施の形態１の画像処理装置１は、重み係数算出処理および重み付き色素量推定処理を実行する。なお、ここで説明する処理は、記憶部１５０に格納された画像処理プログラム１５１に従って画像処理装置１の各部が動作することによって実現される。
【００４３】
先ず、重み係数算出処理について、図６を参照して説明する。図６に示すように、重み係数算出処理では、先ず、推定オペレータ算出部１４１が、推定オペレータＷを算出する（ステップＳ１１）。具体的には、背景技術で示した次式（５）に従って推定オペレータＷを算出する。
【数８】
【００４４】
ここで、背景技術で示したように、次式（３）で定義されるシステム行列Ｈを導入する。
Ｈ＝ＦＳＥ ・・・（３）
光学フィルタ１１９１ａ，１１９１ｂの分光透過率Ｆ、ＲＧＢカメラ１１１の分光感度特性Ｓおよび単位時間当たりの照明の分光放射特性Ｅ（＾）は、使用する機器を選定の後、分光計等を用いて予め測定しておく。なお、ここでは、光学系１１７の分光透過率は１.０と近似しているが、この近似値１.０からの乖離が許容できない場合には、光学系１１７の分光透過率も予め測定し、照明の分光放射特性Ｅに乗じればよい。また、標本の分光透過率の自己相関行列ＲSSおよびＲＧＢカメラ１１１のノイズの自己相関行列ＲNNについても、事前に測定しておく。ＲSSは、Ｈ＆Ｅ染色された典型的な標本を用意し、分光計によって複数の点の分光透過率を測定して自己相関行列を求めることによって得られる。ＲNNは、標本無しの状態で画像取得部１１０によってマルチバンド画像を取得し、得られた６バンドのマルチバンド画像の各バンドについて画素値の分散を求め、これを対角成分とする行列を生成することによって得られる。ただし、バンド間でノイズの相関はないと仮定している。
【００４５】
図７は、横軸を波長λ（ｎｍ）とし、縦軸を波長λにおけるバンドｂの推定オペレータ値ｗλbとして、推定オペレータＷの各波長λに対応する要素の値（推定オペレータ値ｗλb）をバンド毎（バンド１〜バンド６）にプロットして推定オペレータＷをグラフ化した図である。算出された推定オペレータＷのデータは、記憶部１５０に保持される。
【００４６】
次に、図６に示すように、影響度算出部１４３が、ステップＳ１１で算出された推定オペレータＷからバンド毎のノイズに対する相対的な影響度を算出する（ステップＳ１３）。具体的には、次式（１０）に従ってバンド毎のノイズに対する相対的な影響度を算出する。
【数９】
Ｂはバンド数を表す。ｃλbは、推定オペレータ値ｗλbのノイズに対する相対的な影響度である。ｃλbの値は、波長λにおけるバンドｂの推定オペレータ値ｗλbを、波長λにおける全バンドの推定オペレータ値の合計値で割ることによって得た値（波長λにおける全バンドの推定オペレータ値に対する所定のバンドｂにおける推定オペレータ値ｗλbの相対的な大きさ）であり、すなわち、推定オペレータＷを用いて推定される染色標本の分光透過率の波長λにおける推定値に対する推定オペレータ値ｗλbの相対的な影響度を意味する。ここで、推定オペレータ値ｗλbは、入力値が原信号でも観測ノイズでも同様に影響する。したがって、ｃλbの値は、バンド毎の観測ノイズに対する相対的な影響度を意味する。
【００４７】
図８は、横軸を波長λ（ｎｍ）、縦軸をバンド毎のノイズに対する相対的な影響度ｃλbとして、算出されたバンド毎のノイズに対する相対的な影響度ｃλbの一例をグラフ化した図である。図８に示す例では、５４０〜５６０ｎｍの波長帯域におけるバンド毎のノイズに対する相対的な影響度ｃλbの値が、他の波長帯域と比較して高い値を示している。このように、バンド毎のノイズに対する相対的な影響度ｃλbは、ノイズ（観測ノイズ）に対する波長依存性を示しており、このバンド毎のノイズに対する相対的な影響度ｃλbを算出して解析することによって、推定オペレータＷを用いたスペクトル推定における観測ノイズの影響度、すなわち推定オペレータＷのノイズ感度を予測できる。具体的には、その値の大きい波長帯域では観測ノイズによる影響を受け易く、観測ノイズに対するノイズ耐性が低いと予測できる。一方、バンド毎のノイズに対する相対的な影響度ｃλbが小さい波長帯域では観測ノイズによる影響を受け難く、観測ノイズに対するノイズ耐性が高いと予測できる。すなわち、図８の例では、５４０〜５６０ｎｍの波長帯域が、他の波長帯域と比較して観測ノイズによる影響を受け易いと予測できる。実施の形態１では、このバンド毎のノイズに対する相対的な影響度ｃλbを用いて観測ノイズによる影響を受け難い波長帯域におけるスペクトル推定の結果が色素量推定に反映され易く、観測ノイズによる影響を受け易い波長帯域におけるスペクトル推定の結果が色素量推定に反映され難くなるように重み係数ωλを波長毎に算出する。
【００４８】
すなわち、図６に示すように、重み算出部１４５が、対象標本の色素量を推定する際に用いる重み係数を波長毎に算出する（ステップＳ１５）。具体的には先ず、次式（１１）に従って、バンド毎のノイズに対する相対的な影響度ｃλbから波長毎のノイズに対する相対的な影響度ｃλを算出する。
【数１０】
なお、ここでは、先鋭的特徴を捉えるため、バンド毎のノイズに対する相対的な影響度ｃλbの絶対値の最大値を用いたが、絶対値の平均値やノルムを用いて波長毎のノイズに対する相対的な影響度ｃλを求めることとしてもよい。例えば、図８に示す例においてバンド２に着目すれば、５２０〜５４０ｎｍの波長帯域で得られたプロットＰ１１に対応するノイズに対する相対的な影響度ｃλbがノイズに対する相対的な影響度ｃλとして算出される。また、バンド５に着目すれば、５４０〜５６０ｎｍの波長帯域で得られたプロットＰ１３に対応するノイズに対する相対的な影響度ｃλbがノイズに対する相対的な影響度ｃλとして算出される。
【００４９】
次いで、波長λにおける重み係数ωλを算出する。ここで、色素量推定の際に、ノイズに対する相対的な影響度ｃλの低い波長帯域に設定される重み係数ωλが大きく、ノイズに対する相対的な影響度の高い波長帯域に設定される重み係数ωλが小さくなるように波長λにおける重み係数ωλを算出する。例えば、次式（１２），（１３），（１４）に従って波長λにおける重み係数ωλを算出し、波長λにおけるノイズに対する相対的な影響度ｃλが低い場合には重み係数ωλの値を大きくし、ノイズに対する相対的な影響度ｃλが高い場合には重み係数ωλの値を小さくする。
【数１１】
【００５０】
また、αは所定の閾値であり、このαの値によって、ノイズに対する相対的な影響度ｃλの高い波長帯域を色素量の推定にどこまで利用するかを決定する。αの値を高くしすぎると重み係数ωλが常に１となり、重み係数ωλの役割を果たさなくなる。一方、αの値を低くしすぎるとｃλ＞１の波長帯域全てに重み係数ωλが設定されてしまう。例えば、図８に示した例では、ほとんどの波長帯域に重み係数ωλが設定されることとなる。したがって、予めノイズに対する相対的な影響度を考慮せずに色素量推定を行った結果から決定する等して最も推定精度が高まるようにαの値を設定する必要がある。
【００５１】
なお、ここでは、ノイズに対する相対的な影響度ｃλが相対的に高い波長帯域の重み係数ωλに０以上１以下の値を設定することとしたが、０を設定することによって、ノイズに対する相対的な影響度ｃλが相対的に高い波長帯域を色素量の推定に用いないこととしてもよい。また、ノイズに対する相対的な影響度ｃλが相対的に高い波長帯域にのみ重み係数ωλを設定することとしたが、全ての波長帯域に対して相対的に重み係数ωλを設定してもよい。算出された重み係数ω1，ω2，・・・，ωDは、記憶部１５０に保持される。
【００５２】
そして、次式（１５）に従って、推定オペレータＷのノイズに対する相対的な影響度ｃを算出する。
【数１２】
なお、ここでは、先鋭的特徴を捉えるため、ｃ1，ｃ2，・・・，ｃDの最大値を用いたが、平均値やノルムを用いて推定オペレータＷのノイズに対する相対的な影響度ｃを求めることとしてもよい。例えば、図８に示す例では、５４０〜５６０ｎｍの波長帯域で得られたプロットＰ１３に対応するバンド５のノイズに対する相対的な影響度ｃλが、ノイズに対する相対的な影響度ｃとして算出される。
【００５３】
ここで求めたノイズに対する相対的な影響度ｃは、推定オペレータＷの信頼性を表す尺度として用いることができ、例えば画像処理装置１の装置信頼性を評価するために参照される。すなわち、ノイズに対する相対的な影響度ｃの値が小さい程スペクトル推定における観測ノイズの影響度は低く、推定オペレータＷの信頼性が高い。一方、ノイズに対する相対的な影響度ｃの値が大きい場合には、スペクトル推定における観測ノイズの影響度が高くなり、推定オペレータＷの信頼性は低下する。このノイズに対する相対的な影響度ｃを用いて画像処理装置１の装置信頼性を評価する際には、その値が小さければ画像処理装置１の装置信頼性を高いと評価する。一方、ノイズに対する相対的な影響度ｃの値が大きければ画像処理装置１の装置信頼性は低いと評価する。
【００５４】
続いて、重み付き色素量推定処理の手順について、図９を参照して説明する。図９に示すように、重み付き色素量推定処理では、先ず、マルチバンド画像取得制御部１６１が、画像取得部１１０の動作を制御して色素量の推定対象の対象標本をマルチバンド撮像し、対象標本画像を取得する（ステップＳ２１）。
【００５５】
次に、スペクトル推定部１４７が、ステップＳ２１で取得した対象標本画像の画素値をもとに、対象標本のスペクトル（分光透過率）を推定する（ステップＳ２３）。具体的には、図６のステップＳ１１で算出した推定オペレータＷを用いる。そして、背景技術で示した次式（４）に従い、対象標本画像の推定対象画素である任意の点ｘにおける画素の画素値の行列表現Ｇ（ｘ）から、対応する対象標本の標本点における分光透過率の推定値Ｔ＾（ｘ）を推定する。得られた分光透過率の推定値Ｔ＾（ｘ）は、記憶部１５０に格納される。
【数１３】
【００５６】
そして、色素量推定部１４９が、ステップＳ２３で推定した分光透過率の推定値Ｔ＾（ｘ）に基づき、重み係数算出処理（図６）で算出した重み係数ωλを用いて対象標本の色素量を推定する（ステップＳ２５）。ここで、色素量推定部１４９は、対象標本の染色に用いたヘマトキシリンおよびエオジンの各色素の基準分光特性をもとに、対象標本画像の任意の点ｘに対応する標本点における細胞核を染色したヘマトキシリン（色素Ｈ）の色素量、細胞質を染色したエオジン（色素Ｅ）の色素量、および赤血球を染色したエオジン（色素Ｒ）の色素量を推定する。具体的には、対象標本画像の点ｘにおける分光透過率の推定値Ｔ＾（ｘ）に基づいて、点ｘに対応する対象標本の標本点に固定された色素Ｈ，色素Ｅ，色素Ｒそれぞれの色素量を推定する。すなわち、背景技術で示した次式（９）を複数の波長λに関して連立させ、ｄH，ｄE，ｄRについて解く。
【数１４】
【００５７】
例として、３つの波長λ1，λ2，λ3について式（９）を連立させた場合を考えると、次式（１６）のように行列表記できる。
【数１５】
ここで、式（１６）を次式（１７）に置き換える。
【数１６】
この関係式（１７）から、最小二乗法を用いて色素量ｄH，ｄE，ｄRを算出する。最小二乗法とは、単回帰式において誤差の二乗和を最小にするようにｄ（ｘ）を決定する方法であり、一般的には次式（１８）で算出できる。
【数１７】
【００５８】
実施の形態１では、重み係数算出処理で算出した重み係数ωλ（ω1，ω2，・・・，ωD）を加味して色素量ｄH，ｄE，ｄRを算出する。このため、次式（１９），（２０）を用いる。
【数１８】
ωは、重み係数ωλに対応するＤ行Ｄ列の行列であり、ｄｉａｇ（）は、対角行列を表す。以上のようにして重み係数ωλを用いて推定された対象標本画像の点ｘに対応する染色標本の標本点における色素量ｄH，ｄE，ｄRは、記憶部１５０に保持される。そして、例えば、ここで推定された色素量ｄH，ｄE，ｄRに基づいて対象標本画像の色が補正され、カメラの特性や染色状態のばらつき等が補正されて表示用のＲＧＢ画像が合成される。このＲＧＢ画像は、表示部１３０に画面表示されて病理診断に利用される。
【００５９】
実施の形態１によれば、推定オペレータＷからバンド毎のノイズに対する相対的な影響度ｃλbを算出し、波長毎のノイズに対する相対的な影響度ｃλを算出することができ、推定オペレータＷを用いたスペクトル推定における観測ノイズの影響度、すなわち推定オペレータＷのノイズ感度を予測できる。また、波長毎のノイズに対する相対的な影響度ｃλから波長毎に重み係数ωλを算出し、算出した重み係数ωλを用いて対象標本の色素量を推定することができる。具体的には、ノイズに対する相対的な影響度ｃλの低い波長帯域の重み係数ωλを大きく、ノイズに対する相対的な影響度ｃλの高い波長帯域の重み係数ωλを小さくして色素量の推定を行うことができるので、推定オペレータＷを用いたスペクトル推定時に受ける観測ノイズの影響を加味して色素量を推定することができ、推定精度の低下を軽減できる。また、ノイズに対する相対的な影響度ｃλをもとに算出したノイズに対する相対的な影響度ｃを用いて画像処理装置１の装置信頼性を評価することができ、推定オペレータＷを用いたスペクトル推定時に受ける観測ノイズの影響下における装置の信頼性を適正に判断することができる。
【００６０】
なお、実施の形態１では、画像処理装置１において、推定オペレータ算出部１４１が推定オペレータＷを算出し、重み算出部１４５が重み係数ωλを算出することによって重み係数ωλを取得する場合について説明したが、推定オペレータＷや重み係数ωλは、推定の都度毎回算出しなくてもよい。すなわち、予め推定オペレータＷおよび／または重み係数ωλを算出し、記憶部１５０に格納しておく構成としてもよい。そして、スペクトル推定の際には記憶部１５０から推定オペレータＷのデータを読み出して用い、色素量推定の際には記憶部１５０から重み係数ωλのデータを読み出して用いる構成としてもよい。
【００６１】
（実施の形態２）
次に、実施の形態２について説明する。図１０は、実施の形態２の画像処理装置１ｂの機能構成を説明するブロック図である。なお、実施の形態１で説明した構成と同一の構成については、同一の符号を付する。実施の形態２では、画像処理装置１ｂは、図１に示して説明した画像取得部１１０と、入力部１２０と、表示部１３０と、演算部１４０ｂと、記憶部１５０ｂと、装置各部を制御する制御部１６０とを備える。そして、演算部１４０ｂは、影響度算出部１４３ｂと、重み算出部１４５ｂと、スペクトル推定部１４７ｂと、色素量推定部１４９ｂとを含む。
【００６２】
また、記憶部１５０ｂには、予め求められた推定オペレータＷ（ウィナー推定行列）のデータが推定オペレータデータ１５３ｂとして格納されるとともに、推定オペレータデータ１５３ｂに格納された推定オペレータＷと対象標本をマルチバンド撮像した対象標本画像の画素値とからノイズに対する相対的な影響度を解析して色素量推定に用いる重み係数を算出し、この重み係数を用いて対象標本の色素量を推定する処理を実現するための画像処理プログラム１５１ｂが格納される。
【００６３】
図１１は、実施の形態２の画像処理装置１ｂが行う処理手順を示すフローチャートである。なお、ここで説明する処理は、記憶部１５０ｂに格納された画像処理プログラム１５１ｂに従って画像処理装置１ｂの各部が動作することによって実現される。
【００６４】
先ず、マルチバンド画像取得制御部１６１が、画像取得部１１０の動作を制御して色素量の推定対象の対象標本をマルチバンド撮像し、対象標本画像を取得する（ステップＳ３１）。次に、スペクトル推定部１４７ｂが、推定オペレータデータ１５３ｂから推定オペレータＷを読み出し、この推定オペレータＷを用い、ステップＳ３１で取得した対象標本画像の画素値をもとに対象標本のスペクトル（分光透過率）を推定する（ステップＳ３３）。
【００６５】
続いて、影響度算出部１４３ｂが、推定オペレータデータ１５３ｂから読み出した推定オペレータＷと、ステップＳ３１で取得した対象標本画像の画素値とからバンド毎のノイズに対する相対的な影響度を算出する（ステップＳ３５）。具体的には、次式（２１）に従ってバンド毎のノイズに対する相対的な影響度ｃλbを算出する。
【数１９】
ｇ(ｂ)は、バンドｂの対象標本画像の推定対象画素である点ｘの画素値である。このようにして実施の形態２では、推定オペレータ値ｗλbに加えて、推定対象画素の画素値を用いてバンド毎のノイズに対する相対的な影響度ｃλbを算出する。
【００６６】
そして、重み算出部１４５ｂが、波長毎の重み係数を算出する（ステップＳ３７）。重み係数の算出の手法は実施の形態１と同様であり、ステップＳ３５で求めたバンド毎のノイズに対する相対的な影響度ｃλbから式（１１）に従って波長毎のノイズに対する相対的な影響度ｃλを算出し、波長毎のノイズに対する相対的な影響度ｃλから式（１２）〜（１４）に従って波長λにおける重み係数ωλを算出する。また、式（１５）に従って、画像処理装置１ｂの装置信頼性を評価するための値である推定オペレータＷのノイズに対する相対的な影響度ｃを算出する。ここで求めたノイズに対する相対的な影響度ｃは、実施の形態１と同様にして画像処理装置１ｂの装置信頼性を評価するために参照される。
【００６７】
そして、色素量推定部１４９ｂが、ステップＳ３３で推定した分光透過率の推定値Ｔ＾（ｘ）に基づき、ステップＳ３７で算出した重み係数ωλを用いて対象標本の色素量を推定する（ステップＳ３９）。色素量の推定手法は、実施の形態１と同様である。
【００６８】
実施の形態２によれば、実施の形態１と同様の効果を奏するとともに、推定オペレータＷとバンド毎の対象標本画像の画素値ｇ(ｂ)とからバンド毎のノイズに対する相対的な影響度ｃλbを算出し、波長毎のノイズに対する相対的な影響度ｃλを算出することができ、推定オペレータＷを用いたスペクトル推定における観測ノイズの影響、すなわち推定オペレータＷのノイズ感度を対象標本画像の画素値を用いて予測することができる。
【００６９】
なお、実施の形態２では、推定対象画素の画素値を用いてバンド毎のノイズに対する相対的な影響度ｃλbを算出する場合について説明したが、例えば、予め複数の画素値についてバンド毎のノイズに対する相対的な影響度ｃλbを解析し、重み係数ωλを算出して記憶部１５０ｂに格納しておくこととしてもよい。そして、色素量の推定時には、推定対象画素の画素値に応じて最適な重み係数ωλを記憶部１５０ｂから読み出して用いることとしてもよい。例えば、所定幅で画素値をサンプリングし、サンプリングした各画素値についてそれぞれ重み係数ωλを算出しておくこととしてもよい。そして、色素量の推定時には、推定対象画素の画素値と最も値が近い画素値の重み係数ωλを用いて色素量の推定を行う。あるいは、対象標本画像に映る細胞核や細胞質、赤血球、背景といった組織毎に重み係数ωλを算出しておくこととしてもよい。具体的には、組織に応じた代表的な画素値を用いて組織毎の重み係数ωλとして算出しておくこととしてもよい。そして、色素量の推定時には、推定対象画素の画素値から対応する標本点の組織を特定し、特定した組織に応じた重み係数ωλを用いて色素量の推定を行う。
【００７０】
また、実施の形態２では、推定パラーメータＷを予め算出しておく場合について説明したが、推定の都度算出することとしてもよい。
【００７１】
また、上記の実施の形態では、病理標本を撮像したマルチバンド画像から分光透過率のスペクトル特徴値を推定する場合について説明したが、分光特性値として、分光反射率のスペクトル特徴値を推定する場合にも、同様に適用できる。
【図面の簡単な説明】
【００７２】
【図１】画像処理装置の構成を示す図である。
【図２】カラーフィルタの配列例およびＲＧＢ各バンドの画素配列を模式的に示す図である。
【図３−１】一の光学フィルタの分光透過率特性を示す図である。
【図３−２】他の光学フィルタの分光透過率特性を示す図である。
【図４】Ｒ，Ｇ，Ｂ各バンドの分光感度の例を示す図である。
【図５】実施の形態１の画像処理装置の機能構成を示すブロック図である。
【図６】重み係数算出処理の処理手順を示すフローチャートである。
【図７】推定オペレータをグラフ化した図である。
【図８】バンド毎のノイズに対する相対的な影響度の一例をグラフ化した図である。
【図９】重み付き色素量推定処理の処理手順を示すフローチャートである。
【図１０】実施の形態２の画像処理装置の機能構成を示すブロック図である。
【図１１】実施の形態２の画像処理装置が行う処理手順を示すフローチャートである。
【図１２】ＲＧＢ画像の一例を示す図である。
【符号の説明】
【００７３】
１ 画像処理装置
１１０ 画像取得部
１２０ 入力部
１３０ 表示部
１４０ 演算部
１４１ 推定オペレータ算出部
１４３ 影響度算出部
１４５ 重み算出部
１４７ スペクトル推定部
１４９ 色素量推定部
１５０ 記憶部
１５１ 画像処理プログラム
１６０ 制御部
１６１ マルチバンド画像取得制御部
【特許請求の範囲】
【請求項１】
推定オペレータを用いて分光特性の推定を行う画像処理装置において、
前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析手段を備えることを特徴とする画像処理装置。
【請求項２】
染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定手段と、
分光特性の波長毎の重み係数を取得する重み係数取得手段と、
前記分光特性推定手段によって推定された分光特性と前記重み係数取得手段によって取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定手段と、
を備えることを特徴とする画像処理装置。
【請求項３】
前記重み係数取得手段は、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析することで、前記波長毎の重み係数を算出することを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記影響度解析手段によって解析されたノイズに対する相対的な影響度から波長毎の重み係数を設定する重み係数設定手段を備えることを特徴とする請求項１に記載の画像処理装置。
【請求項５】
染色標本を撮像した染色標本画像をもとに、前記推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定手段と、
前記分光特性推定手段によって推定された分光特性と、前記重み係数設定手段によって設定された重み係数とをもとに、前記染色標本の色素量を推定する色素量推定手段と、
を備えることを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記重み係数設定手段は、前記ノイズに対する相対的な影響度の高い波長に対して設定する重み係数を小さくすることを特徴とする請求項４または５に記載の画像処理装置。
【請求項７】
前記重み係数設定手段は、前記ノイズに対する相対的な影響度の高い波長に対して設定する重み係数を０にすることを特徴とする請求項４〜６のいずれか１つに記載の画像処理装置。
【請求項８】
前記影響度解析手段は、前記染色標本画像の画素値を用いて前記ノイズに対する相対的な影響度を解析することを特徴とする請求項１に記載の画像処理装置。
【請求項９】
推定オペレータを用いて分光特性の推定を行う画像処理装置の信頼性評価方法であって、
前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析し、解析したノイズに対する相対的な影響度をもとに前記画像処理装置の装置信頼性を評価することを特徴とする画像処理装置の信頼性評価方法。
【請求項１０】
推定オペレータを用いて分光特性の推定を行うコンピュータに、
前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析ステップを実行させることを特徴とする画像処理プログラム。
【請求項１１】
コンピュータに、
染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定ステップと、
分光特性の波長毎の重み係数を取得する重み係数取得ステップと、
前記分光特性推定ステップで推定された分光特性と前記重み係数取得ステップで取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定ステップと、
を実行させることを特徴とする画像処理プログラム。
【請求項１】
推定オペレータを用いて分光特性の推定を行う画像処理装置において、
前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析手段を備えることを特徴とする画像処理装置。
【請求項２】
染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定手段と、
分光特性の波長毎の重み係数を取得する重み係数取得手段と、
前記分光特性推定手段によって推定された分光特性と前記重み係数取得手段によって取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定手段と、
を備えることを特徴とする画像処理装置。
【請求項３】
前記重み係数取得手段は、前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析することで、前記波長毎の重み係数を算出することを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記影響度解析手段によって解析されたノイズに対する相対的な影響度から波長毎の重み係数を設定する重み係数設定手段を備えることを特徴とする請求項１に記載の画像処理装置。
【請求項５】
染色標本を撮像した染色標本画像をもとに、前記推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定手段と、
前記分光特性推定手段によって推定された分光特性と、前記重み係数設定手段によって設定された重み係数とをもとに、前記染色標本の色素量を推定する色素量推定手段と、
を備えることを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記重み係数設定手段は、前記ノイズに対する相対的な影響度の高い波長に対して設定する重み係数を小さくすることを特徴とする請求項４または５に記載の画像処理装置。
【請求項７】
前記重み係数設定手段は、前記ノイズに対する相対的な影響度の高い波長に対して設定する重み係数を０にすることを特徴とする請求項４〜６のいずれか１つに記載の画像処理装置。
【請求項８】
前記影響度解析手段は、前記染色標本画像の画素値を用いて前記ノイズに対する相対的な影響度を解析することを特徴とする請求項１に記載の画像処理装置。
【請求項９】
推定オペレータを用いて分光特性の推定を行う画像処理装置の信頼性評価方法であって、
前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析し、解析したノイズに対する相対的な影響度をもとに前記画像処理装置の装置信頼性を評価することを特徴とする画像処理装置の信頼性評価方法。
【請求項１０】
推定オペレータを用いて分光特性の推定を行うコンピュータに、
前記推定オペレータから分光特性の波長軸におけるノイズに対する相対的な影響度を解析する影響度解析ステップを実行させることを特徴とする画像処理プログラム。
【請求項１１】
コンピュータに、
染色標本を撮像した染色標本画像をもとに、推定オペレータを用いて前記染色標本の分光特性を推定する分光特性推定ステップと、
分光特性の波長毎の重み係数を取得する重み係数取得ステップと、
前記分光特性推定ステップで推定された分光特性と前記重み係数取得ステップで取得された波長毎の重み係数とをもとに、前記染色標本の色素量を推定する色素量推定ステップと、
を実行させることを特徴とする画像処理プログラム。
【図１】
【図２】
【図３−１】
【図３−２】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１２】
【図１１】
【図２】
【図３−１】
【図３−２】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１２】
【図１１】
【公開番号】特開２００９−１８６４４２（Ｐ２００９−１８６４４２Ａ）
【公開日】平成２１年８月２０日（２００９．８．２０）
【国際特許分類】
【出願番号】特願２００８−２９６９３（Ｐ２００８−２９６９３）
【出願日】平成２０年２月８日（２００８．２．８）
【出願人】（００００００３７６）オリンパス株式会社 (11,466)
【Ｆターム（参考）】
【公開日】平成２１年８月２０日（２００９．８．２０）
【国際特許分類】
【出願日】平成２０年２月８日（２００８．２．８）
【出願人】（００００００３７６）オリンパス株式会社 (11,466)
【Ｆターム（参考）】
[ Back to top ]