JP6215894B2 - 画像処理方法およびシェーディング基準データ作成方法 - Google Patents

Description

この発明は、撮像により得られた画像における明暗ムラを解消するために画像を補正する技術に関するものである。

照明下で撮像対象物を撮像する技術においては、撮像対象物の背景が略一様な光学濃度を有するものであっても、例えば撮像対象物に入射する照明光の光量ムラなど、照明系または撮像系の特性に起因する明暗ムラが画像に現れることがある。このような明暗ムラを補正するための技術としてシェーディング補正技術がある。例えば特許文献１に記載の技術では、顕微鏡を用いて標本（撮像対象物）を撮像した画像の各画素に対し、標本なしで撮像された画像のシェーディング特性に基づいて作成されたシェーディング補正データを作用させることにより、光学系に起因する画像濃度のムラの解消が図られている。

特開２００５−０７２９６７号公報

例えば液体中の撮像対象物を撮像する場合、撮像対象物に対し照明光が液面を介して入射することになるが、液面の状態は必ずしも一定ではない。この事例のように、撮像対象物ごとに照明光の入射条件が異なる場合がある。このような場合、事前に取得したシェーディング特性が必ずしも適切なものとならず、有効な補正を行えないことがあり得る。

この問題に関して上記した従来技術では、モニタに表示される補正後の画像を確認しながら操作者が補正特性の調整を行うことができるように構成されている。しかしながら、調整作業を効率的に行うには専門的な知識を必要とし、また補正結果が操作者の主観に左右されるため、補正結果の安定性という点で問題がある。

この発明は上記課題に鑑みなされたものであり、撮像により得られた画像における明暗ムラを効果的に、かつ安定的に補正することのできる技術を提供することを目的とする。

この発明の一の態様は、上記目的を達成するため、撮像対象物が略一様な背景と共に撮像された原画像を取得する画像取得工程と、前記原画像の二次元輝度プロファイルを二次元正規分布の確率密度関数に近似した近似式を特定する近似工程と、前記原画像を構成する各画素の輝度値を、前記近似式により表される二次元輝度プロファイルの輝度値に基づいて補正する補正工程とを備えた画像処理方法である。

また、この発明の他の態様は、撮像対象物が略一様な背景と共に撮像された原画像に対しシェーディング補正を行う際の基準となるシェーディング基準データの作成方法であって、上記目的を達成するため、前記原画像を取得する画像取得工程と、前記原画像の二次元輝度プロファイルを二次元正規分布の確率密度関数に近似した近似式を特定する近似工程と、前記近似式により表される二次元輝度プロファイルを有する仮想的な画像に対応するデータを前記シェーディング基準データとして設定する設定工程とを備えている。

このように構成された発明では、撮像対象物の背景部分が有する輝度プロファイルが二次元正規分布に従うものと見なされ、該分布を表す確率密度関数の近似式が原画像の二次元輝度プロファイルに基づき特定される。そして、この近似式で表される二次元輝度プロファイルが、シェーディング特性を近似的に表す数学モデルとしての意味を持つ。本願発明者らの知見によれば、このような近似は、例えば液体中の撮像対象物のように、背景が略一様となった多くの撮像において有効である。すなわち、従来技術では事前撮像により取得される光学系のシェーディング特性に対応するプロファイルに代えて、上記のようにして得られた近似式が表す輝度プロファイルを用いて、シェーディング補正技術と同様の補正が行われる。これにより、原画像における明暗ムラを解消して明るさが均一化された画像を得ることができる。そして、処理は操作者の主観的判断に左右されないため、安定した結果を得ることが可能である。

上記のように、本発明の画像処理方法によれば、撮像により得られた画像における明暗ムラを効果的に、かつ安定的に補正することができる。また、本発明のシェーディング基準データの作成方法によれば、画像における明暗ムラを補正するのに好適なシェーディング基準データを作成することができる。

本発明にかかる画像処理方法を適用可能な撮像装置の概略構成を示す図である。 ウェル画像の輝度プロファイルの例を示す図である。 二次元正規分布の確率密度関数を示す図である。 シェーディング補正処理を示すフローチャートである。 モデル化処理の例を示すフローチャートである。 処理対象画像から漸化式の初期値を求める方法を説明する図である。 二値化画像の例を示す図である。 補正画像の例を示す図である。

図１は本発明にかかる画像処理方法を適用可能な撮像装置の概略構成を示す図である。この撮像装置１は、ウェルプレートＷＰの上面に形成されたウェルＷと称される凹部に注入された液体中で培養される細胞、細胞コロニー、細菌等（以下、「細胞等」と称する）の生試料を撮像する装置である。

ウェルプレートＷＰは、創薬や生物科学の分野において一般的に使用されているものであり、平板状のプレートの上面に、断面が略円形の筒状に形成され底面が透明で平坦なウェルＷが複数設けられている。１つのウェルプレートＷＰにおけるウェルＷの数は任意であるが、例えば９６個（１２×８のマトリクス配列）のものを用いることができる。各ウェルＷの直径および深さは代表的には数ｍｍ程度である。なお、この撮像装置１が対象とするウェルプレートのサイズやウェルの数はこれらに限定されるものではなく任意であり、例えば６ないし３８４穴のものが一般的に使用されている。また、複数ウェルを有するウェルプレートに限らず、例えばディッシュと呼ばれる平型の容器で培養された細胞等の撮像にも、この撮像装置１を使用することが可能である。

ウェルプレートＷＰの各ウェルＷには、培地Ｍとしての液体が所定量注入され、この液体中において所定の培養条件で培養された細胞等が、この撮像装置１の撮像対象物となる。培地は適宜の試薬が添加されたものでもよく、また液状でウェルＷに投入された後ゲル化するものであってもよい。後述するように、この撮像装置１では、例えばウェルＷの内底面で培養された細胞等を撮像対象とすることができる。常用される一般的な液量は、５０ないし２００マイクロリットル程度である。

撮像装置１は、ウェルプレートＷＰを保持するホルダ１１と、ホルダ１１の上方に配置される照明部１２と、ホルダ１１の下方に配置される撮像部１３と、これら各部の動作を制御するＣＰＵ１４１を有する制御部１４とを備えている。ホルダ１１は、試料を培地Ｍとともに各ウェルＷに担持するウェルプレートＷＰの下面周縁部に当接してウェルプレートＷＰを略水平姿勢に保持する。

照明部１２は、ホルダ１１により保持されたウェルプレートＷＰに向けて照明光を出射する。照明光の光源としては、例えば白色ＬＥＤ（Light Emitting Diode）を用いることができる。光源と適宜の照明光学系とを組み合わせたものが、照明部１２として用いられる。照明部１２により、ウェルプレートＷＰに設けられたウェルＷ内の細胞等が上方から照明される。

ホルダ１１により保持されたウェルプレートＷＰの下方に、撮像部１３が設けられる。撮像部１３には、ウェルプレートＷＰの直下位置に撮像光学系が配置されており、撮像光学系の光軸は鉛直方向に向けられている。図１は側面図であり、図の左右方向が水平方向を、また上下方向が鉛直方向をそれぞれ表す。

撮像部１３により、ウェルＷ内の細胞等が撮像される。具体的には、照明部１２から出射されウェルＷの上方から液体に入射した光が撮像対象物を照明し、ウェルＷ底面から下方へ透過した光が、撮像部１３の対物レンズ１３１を含む撮像光学系を介して撮像素子１３２の受光面に入射する。撮像光学系により撮像素子１３２の受光面に結像する撮像対象物の像が、撮像素子１３２により撮像される。撮像素子１３２は、一次元の受光面を有するリニアイメージセンサまたは二次元の受光面を有するエリアイメージセンサであり、例えばＣＣＤセンサまたはＣＭＯＳセンサを用いることができる。

撮像部１３は、制御部１４に設けられたメカ制御部１４６により水平方向および鉛直方向に移動可能となっている。具体的には、メカ制御部１４６がＣＰＵ１４１からの制御指令に基づき駆動機構１５を作動させ、撮像部１３を水平方向に移動させることにより、撮像部１３がウェルＷに対し水平方向に移動する。また鉛直方向への移動によりフォーカス調整がなされる。

また、駆動機構１５は、撮像部１３を水平方向に移動させる際、図において点線矢印で示すように照明部１２を撮像部１３と一体的に移動させる。すなわち、照明部１２は、その光中心が撮像部１３の光軸と略一致するように配置されており、撮像部１３が水平方向に移動するとき、これと連動して移動する。

１つのウェルＷについて詳細な画像が必要とされるとき、撮像視野内に１つのウェルＷの全体が収められた状態で撮像が行われる。このとき、メカ制御部１４６は、光軸が当該ウェルＷの中心と一致するように、撮像部１３を水平方向に位置決めする。これにより、どのウェルＷが撮像される場合でも、当該ウェルＷの中心および照明部１２の光中心が常に撮像部１３の光軸上に位置することとなり、各ウェルＷに対する照明条件を一定にして、撮像条件を良好に維持することができる。

一方、複数のウェルＷについてより高速に撮像することが必要とされるときには、水平方向に撮像部１３を走査移動させながら撮像が行われる。すなわち、メカ制御部１４６が照明部１２および撮像部１３を一体的に一定速度で水平移動させることで走査移動が実現され、撮像部１３が定期的に撮像を行うことで広範囲が撮像される。特に撮像素子１３２としてのリニアイメージセンサを有する構成においては、このような撮像対象物に対する撮像部１３の走査移動により、撮像対象物の二次元画像が取得される。

撮像部１３の撮像素子１３２から出力される画像信号は、制御部１４に送られる。すなわち、画像信号は制御部１４に設けられたＡＤコンバータ（Ａ／Ｄ）１４３に入力されてデジタル画像データに変換される。ＣＰＵ１４１は、受信した画像データに基づき適宜画像処理を実行する。制御部１４はさらに、画像データを記憶保存するための画像メモリ１４４と、ＣＰＵ１４１が実行すべきプログラムやＣＰＵ１４１により生成されるデータを記憶保存するためのメモリ１４５とを有しているが、これらは一体のものであってもよい。ＣＰＵ１４１は、メモリ１４５に記憶された制御プログラムを実行することにより、後述する各種の演算処理を行う。

その他に、制御部１４には、インターフェース（ＩＦ）部１４２が設けられている。インターフェース部１４２は、ユーザからの操作入力の受け付けや、ユーザへの処理結果等の情報提示を行うユーザインターフェース機能のほか、通信回線を介して接続された外部装置との間でのデータ交換を行う機能を有する。ユーザインターフェース機能を実現するために、インターフェース部１４２には、ユーザからの操作入力を受け付ける入力受付部１４７と、ユーザへのメッセージや処理結果などを表示出力する表示部１４８とが接続されている。

なお、制御部１４は、上記したハードウェアを備えた専用装置であってもよく、またパーソナルコンピュータやワークステーション等の汎用処理装置に、後述する処理機能を実現するための制御プログラムを組み込んだものであってもよい。すなわち、この撮像装置１の制御部１４として、汎用のコンピュータ装置を利用することが可能である。汎用処理装置を用いる場合、撮像装置１には、撮像部１３等の各部を動作させるために必要最小限の制御機能が備わっていれば足りる。

図２はウェル画像の輝度プロファイルの例を示す図である。図２上段に示すように、撮像されるウェルＷには培地Ｍが注入されており、培地Ｍの内部、一般的にはウェルＷの内底面Ｗｂに撮像対象物である細胞等Ｃが分布している。培地Ｍは液体、または液体として注入された後ゲル化したものであり、その液面は必ずしも水平面とならず、一般的には下に凸のメニスカスを形成している。

このため、上方から液面に入射する照明光Ｌｉは、ウェルＷの中心に近い部分ではほぼ直進し、ウェル内底部Ｗｂを介してウェルＷの下方へ透過するのに対して、ウェルＷの側壁面に近い位置では液面で屈折し、その進行方向がウェルＷの中心から外方向へ向かう向きに変えられる。その結果、撮像された画像では、図２下段の輝度プロファイル例に示すように、ウェルＷ中心付近では比較的高輝度が得られる一方、ウェルＷ周縁部に近づくにつれて輝度が低下する。この他、培地Ｍ中の細胞等Ｃや画像のノイズ等に起因する輝度の小刻みな変動が画像に現れ得る。

培地Ｍそのものは略均一であり、その光学濃度も概ね一様であると考えられる。したがって細胞等Ｃの背景となる培地Ｍの画像濃度も、本来的には略均一となるはずである。しかしながら、上記した理由により、実際の画像では例えば中央部が明るく周縁部で暗いという明暗ムラが生じる。また、照明光の入射方向の傾きによっても明暗ムラが生じ得る。

このような明暗ムラを補正する技術がシェーディング補正技術である。しかしながら、液体注入時のばらつき等により培地Ｍの量や表面状態がウェルＷごとに必ずしも一定とならず、また照明光の入射光量分布がウェルＷごとに異なることもある。これらの原因により、一般的なシェーディング補正技術を適用することが適当でない場合があり得る。すなわち、従来のシェーディング補正技術においては基準となる照明条件で取得された基準画像が必要であるが、上記のようなケースでは一義的に基準を定めることができない。

そこで、この実施形態では、撮像対象物を内包するウェルＷを撮像して得られた画像を原画像として、原画像から基準画像に対応するデータを抽出し補正を行う。具体的には、培地Ｍの表面状態や照明光量分布に起因する明暗ムラが二次元正規分布に従う輝度分布を有すると仮定して、この輝度分布を近似的に表す確率密度関数を特定する。そして、求められた確率密度関数で表される輝度分布を有する仮想的なモデル画像を基準画像に代わるものとして原画像を補正することで、原画像における明暗ムラが解消された画像を得る。

図３は二次元正規分布の確率密度関数を示す図である。より詳しくは、図３（ａ）は二次元正規分布に従う輝度プロファイルを三次元マップとして表したものであり、図３（ｂ）はこれを等輝度線により二次元マップとして表したものである。画像の横方向座標をＸ座標、縦方向座標をＹ座標としたとき、画像を構成する各画素のうちＸ座標値がｘ、Ｙ座標値がｙで表される画素Ｐ（ｘ，ｙ）の輝度値Ｌ（ｘ，ｙ）については、二次元正規分布に対応する確率密度関数である下記の（式１）により表すことができる。

上記（式１）におけるパラメータａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇの値により輝度値Ｌ（ｘ，ｙ）は種々の値を取り得るが、図３（ａ）および図３（ｂ）のグラフ形状から明らかなように、ＸＹ平面内の任意の直線上での輝度プロファイルは正規分布曲線となる。また、等輝度面は一般に楕円形となり、その長軸の方向も様々な方向を取り得る。

（式１）と図３（ａ）に示される曲面とを対比すると、パラメータａは曲面のピークにおける輝度値に対応し、パラメータｂ，ｃは曲面の広がり形状に対応する。また、パラメータｅ，ｆはピークの座標位置に対応し、パラメータｄは楕円主軸の回転を表す。また、パラメータｇは曲面の裾が漸近する輝度値に対応する。

実データの分布を近似的に表す数式を求める方法としては最小二乗法が広く用いられるが、単に誤差を最小化することのみを目的として、一般的な線形高次多項式による近似表現を許容した場合、ウェル内の輝度分布を表すのに適さない場合が生じ得る。より具体的には、ウェルＷ内を撮像した画像の背景部分においては、ウェルＷの中央付近が最も輝度値が高く、周端部において輝度値は一定値に収束する。このような中央部および終端部の輝度値が拘束条件として近似式に反映されるようになっていなければ、補正後の画像が不自然なものとなってしまうおそれがある。

この点において、近似式を（式１）により表される二次元正規分布の確率密度関数とすることには妥当性がある。また、本願発明者らは、後述するこの実施形態における補正方法を各種の試料を撮像した画像を用いて評価し、多くの場合で良好な結果が得られることを確認している。

撮像条件に起因して原画像に現れる明暗ムラをこのような輝度プロファイルに近似し、この輝度プロファイルの寄与分を原画像から取り除くような補正処理を行うことで、撮像対象物が本来有する光学濃度が反映された画像を得ることができる。このような補正処理は、事前に取得された基準画像に基づく従来技術のシェーディング補正処理に対応するものである。より具体的には、従来技術のシェーディング補正処理において用いられる基準画像に代えて、原画像に基づいて作成されたモデル画像を用いた補正処理である。以下では便宜的に、本実施形態の補正処理をこの実施形態における「シェーディング補正処理」と称することとする。

図４はこの実施形態におけるシェーディング補正処理を示すフローチャートである。この処理は、制御部１４のＣＰＵ１４１がメモリ１４５に予め記憶されている制御プログラムを実行することにより実現される。最初に、ウェルプレートＷＰに設けられたウェルＷをその内容物とともに撮像した原画像が取得される（ステップＳ１０１）。原画像としては撮像部１３によって撮像されたウェルＷの画像を用いることができるほか、外部の撮像装置やデータ記憶装置からインターフェース１４２を介して与えられる画像を用いることができる。取得した画像が複数のウェルＷを含む場合、あるいはウェルＷ以外の余白部分を多く含む場合には、当該画像のうち１つのウェルＷ全体と最小限の余白部分とを含む矩形領域がウェル領域として切り出され、切り出された部分画像が原画像とされる（ステップＳ１０２）。

このようにして得られた原画像に、例えば培地Ｍ中の気泡や異物など撮像対象物とは明らかに異なる不要なオブジェクトが含まれていれば、適宜の画像処理によりこれが除去されて、シェーディング補正処理の対象となる処理対象画像が作成される（ステップＳ１０３）。この場合の画像処理には、例えば画像の「ゴミ取り」と称される公知の各種画像処理アルゴリズムを適用することができる。また、ユーザにより指定されたオブジェクトを除去するような処理が加えられてもよい。このような処理を行うことでシェーディング補正処理後の画像品質をより良好なものとすることができるが、画像がそのような不要オブジェクトを含まない場合も多く、この処理は省かれてもよい。

また、原画像には撮像対象物である細胞等も含まれている。そのうち、ウェル内で特に大きな面積を占めるものについては不要オブジェクトと同様に除去されてもよい。ただし撮像対象物が一般的な細胞であれば培地Ｍとの濃度差が小さいため、そのまま残存させておいても処理結果に及ぼす影響は小さいと考えられる。というのは、以下に説明するように、この実施形態ではウェル内の背景部分の輝度分布を二次元正規分布と見なして処理するため、細胞等に起因する小さな濃度変化は実質的に無視されることになるからである。

次に、作成された処理対象画像に基づき、従来のシェーディング補正処理における基準画像に相当するモデル画像を作成するためのモデル化処理が実行される（ステップＳ１０４）。以下、モデル化処理の原理について説明する。モデル化処理は、処理対象画像中のウェルＷに対応するウェル領域が、所定の輝度分布を有する背景中に撮像対象物たる細胞等が分布したものであると見なし、その背景部分の輝度分布を表す数学モデル（以下、「背景モデル」という）を特定するための処理である。前述したように、この実施形態では、背景部分の輝度分布が（式１）により表される二次元正規分布に従うと仮定している。したがって、背景モデルを特定するためには、処理対象画像における輝度分布によくフィッティングするように、（式１）における未知パラメータａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇを決定すればよい。

フィッティングは、処理対象画像の各画素Ｐ（ｘ，ｙ）の輝度値Ｌ（ｘ，ｙ）と、背景モデルにおいて対応する画素の輝度値Ｆ（ｘ，ｙ）との二乗誤差が最小となるように、未知パラメータａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇを定めることにより行われる。二乗誤差Ｅは下記の（式２）により表すことができる。

（式２）において、符号ｉは画像を構成する個々の画素を区別するための添え字である。処理対象画像はｎ個の画素により構成され、各画素が添え字ｉ（ｉ＝１，２，…，ｎ）により区別される。第ｉの画素Ｐｉの座標値が（ｘ_ｉ，ｙ_ｉ）により表される。また、符号ａ₀〜ｇ₀はパラメータａ〜ｇ各々の初期値としての具体的な数値であることを示す。また符号Δａ〜Δｇはパラメータａ〜ｇ各々の差分を表す。

また、（式２）において符号ΔＦで表される式の各項に含まれる偏微分係数は、（式１）からそれぞれ次のように表すことができる。

（式２）の右辺をΔａ〜Δｇで偏微分するとゼロとなる条件から、次の方程式が得られる。

（式４）において、角括弧を用いて表記される演算子[Ｑ]は、画素ごとに求められる関数Ｑの値を全画素について合計した値を表すものとする。また、符号Δａ₀〜Δｇ₀はパラメータａ〜ｇ各々の初期値からの差分値を示す。また右辺において、
Ｆ_ｉ＝Ｆ（ｘ_ｉ，ｙ_ｉ）
である。（式４）中の偏微分係数の値は、初期値ａ₀〜ｇ₀と各点の座標値（ｘ_ｉ，ｙ_ｉ）とを（式３）に代入することで数値的に求められる。なお、（式４）左辺をなす７行７列の行列の各成分において、一次微分値同士の乗算の順序は入れ換え可能であり、したがって当該行列は対称行列である。

（式４）で示される方程式を解き、初期値ａ₀〜ｇ₀と差分値Δａ₀〜Δｇ₀とを用いてパラメータａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇをそれぞれ求めることで、（式１）により示される背景モデルが確定する。この実施形態では、前記した拘束条件を反映させて特にウェル周端部での背景モデルの輝度分布を実測（すなわち処理対象画像における背景輝度分布）とよく一致させるために、方程式の解法としてＬｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔ法（以下、「ＬＭ法」と略称する）が用いられる。ＬＭ法は、未知パラメータに適宜の初期値を与えた漸化式について反復計算を行って、パラメータを最適値に収束させる計算手法の一例である。拘束条件を反映させた初期値を設定することで、収束結果と実測結果とをよく合致させることができる。

漸化式は、項番号を表す添え字Ｋを用いて、下記の（式５）により表される。

（式５）において、符号Ｈｕは（式４）左辺の７行７列の行列を、また符号∇ｕJは（式４）右辺のベクトルを表すものとする。また、符号Ｄ［Ｈｕ］は行列Ｈｕの対角項を取る関数を表し、符号ｈはＭａｒｑｕａｒｄｔパラメータとも呼ばれる係数である。以下、（式５）を用いて処理対象画像の背景モデルを特定するモデル化処理の具体的手順について説明する。

図５はモデル化処理の例を示すフローチャートである。この処理は、制御部１４のＣＰＵ１４１がメモリ１４５に予め記憶されている制御プログラムを実行することにより実現される。最初に、背景モデルをどのような数学モデルを用いて表すかが確定される（ステップＳ２０１）。一般的には背景モデルとして種々のものが考えられるが、上記の通り、本実施形態では（式１）、すなわち二次元正規分布の確率密度関数が背景モデルとされる。言い換えれば、図５に示す処理は、背景モデルが他の数式で表される場合にも利用可能である。

そして、処理対象画像と背景モデルとの二乗誤差を表す式における各偏微分係数が、背景モデルの式を用いて求められる（ステップＳ２０２）。この実施形態では、（式３）により各偏微分係数が表される。これらの偏微分係数が（式５）に代入される。

次に、（式５）の漸化式に与えられるパラメータの初期値が設定される（ステップＳ２０３）。パラメータａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇの初期値は、処理対象画像を構成する各画素の輝度値に基づき、以下のようにして定められる。

図６は処理対象画像から漸化式の初期値を求める方法を説明する図である。より具体的には、図６（ａ）は処理対象画像の一例を示し、図６（ｂ）は当該画像を構成する画素の輝度値の分布を表すヒストグラムである。また、図６（ｃ）はヒストグラムから求められる各種データを示す図である。図６（ａ）に示すように、中央部が明るく周縁部が暗く写ったウェルＷの画像を処理対象画像Ｉｐとする場合を考える。ここでは細胞等を含まない培地Ｍが注入されたウェルＷを撮像した画像Ｉｐが用いられている。ウェルＷの左上部分に見られる大きな欠損部は液面に形成された気泡に起因するものである。

処理対象画像ＩｐのうちウェルＷに対応する円形部分を構成する各画素における輝度値の分布をヒストグラムで表すと、図６（ｂ）に示すように、比較的低い輝度においてウェルＷ周縁部の低輝度部分に対応するピークが現れ、ウェルＷ内各画素の輝度値はこれより高輝度側に広く分布している。このことは、ウェルＷ内で明暗のムラが大きいことを示している。このヒストグラムから、図６（ｃ）に示す各データが得られる。このうちモードの値２９が、図３（ａ）の曲面の裾部分が漸近する輝度値に対応するパラメータｇの初期値とされる。また、最大値からモードを減じた値１３８が、図３（ａ）の曲面のピーク高さに対応するパラメータａの初期値とされる。

他のパラメータは、図３（ａ）の曲面をＸＹ平面に平行な平面で切断したときの断面が表す楕円の形状に関連するものである。これらの初期値を求めるために、処理対象画像Ｉｐを適宜の閾値で二値化し楕円領域を抽出する。ここでは閾値として、図３（ａ）の曲面が表す二次元正規分布における半値、すなわち先に設定されたパラメータａ，ｇの初期値の中点に相当する値（＝（ａ−ｇ）／２）が用いられる。

図７は二値化画像の例を示す図である。より具体的には、図７（ａ）は図６（ａ）の画像から求められた二値化画像Ｉｂの例を示し、図７（ｂ）はこの二値化画像Ｉｂから得られる各種データを示す図である。図７（ａ）に示すように、処理対象画像Ｉｐから作成される二値化画像Ｉｂでは、処理対象画像Ｉｐのうち閾値以上の輝度値を有する領域が抽出される。このため、主として中央部の高輝度部分が抽出され、こうして抽出された領域の形状が上記した楕円領域の概略形状を表す。したがって、抽出された高輝度領域の形状から残りのパラメータの初期値が求められる。

なお、画像に現れる明暗ムラに起因して、二値化により抽出される高輝度領域は必ずしも単一とならない。そこで、抽出された高輝度領域が複数ある場合には、そのうち面積最大のものに着目することとする。また、例えば図７（ａ）左上に示される気泡部分のように、複数の高輝度領域が全体として１つの連続した領域をなすことが外観上明らかである場合には、それらを連結した１つの高輝度領域と見なすようにしてもよい。

図７（ａ）の点線は、二値化により抽出された高輝度領域Ｒｈに外接する矩形（外接矩形）Ｒを表しており、図７（ｂ）には当該外接矩形Ｒの頂点座標位置およびサイズが示されている。外接矩形Ｒの幅、すなわちＸ方向長さが４５０であり、これがパラメータｂの初期値とされる。一方、外接矩形Ｒの高さ、すなわちＹ方向長さが４６２であり、これがパラメータｃの初期値とされる。また、外接矩形Ｒの重心位置のＸ座標値が３０１、Ｙ座標値が３３５となり、これらがそれぞれパラメータｅ、ｆの初期値とされる。楕円領域の回転に関連するパラメータｄの初期値は０とされる。

各パラメータの初期値を整理すると、
ａ：１３８（ヒストグラムにおける最大値とモードとの差）
ｂ：４５０（二値化画像における外接矩形の幅）
ｃ：４６２（二値化画像における外接矩形の高さ）
ｄ：０（回転なしと仮定）
ｅ：３０１（二値化画像における外接矩形の重心のＸ座標）
ｆ：３３５（二値化画像における外接矩形の重心のＹ座標）
ｇ：２９（モード）
となる。

これらを初期値とするパラメータのセット（ａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇ）を用い、係数ｈを適宜の値、例えば０．０００１に設定して（ステップＳ２０４）、ＬＭ法による演算を実行して（式５）を解くことにより（ステップＳ２０５）、差分値Δａ₀、Δｂ₀、Δｃ₀、Δｄ₀、Δｅ₀、Δｆ₀、Δｇ₀をそれぞれ求める。

求められた差分値Δａ₀、Δｂ₀、Δｃ₀、Δｄ₀、Δｅ₀、Δｆ₀、Δｇ₀をそれぞれパラメータａ、ｂ、ｃ、ｄ、ｅ、ｆ、ｇの現在の設定値から加算または減算して、新たなパラメータセット（ａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇ）を作成する（ステップＳ２０６）。ＬＭ法による反復計算の実行回数が所定値に達していれば（ステップＳ２０７においてＹＥＳ）、処理を終了する。

実行回数が所定値に達していなければ（ステップＳ２０７においてＮＯ）、現在の漸化式の値ｕ^(K)と、新たに計算される値ｕ^(K+1)との差が残差として計算される（ステップＳ２０８）。残差が前回の計算で求められた値より小さければ式が収束に近づいていると判断され（ステップＳ２０９においてＹＥＳ）、１回の計算でパラメータが大きく変動するのを避けるために、係数ｈを現在の値の１０分の１に減少させる（ステップＳ２１０）。一方、残差が前回の計算で求められた値より大きければ式が収束から遠ざかっていると判断され（ステップＳ２０９においてＮＯ）、係数ｈを現在の値の１０倍に増加させる（ステップＳ２１１）。

このようにＬＭ法を用いた反復計算を行うことにより、パラメータセット（ａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇ）が最適化される。その結果、これらのパラメータが代入された（式１）により表される背景モデルが、実際に取得された処理対象画像における輝度分布に応じて最適化される。上記したモデル化処理により作成された背景モデルを仮想的な基準画像と見なし、処理対象画像を構成する各画素の輝度値をモデル画像の画素の輝度値を用いて補正することにより、シェーディング補正処理を施された補正画像が得られる。

図４に戻ってシェーディング補正処理の説明を続ける。モデル化処理が終了すると、処理対象画像の各画素の輝度値が背景モデルに基づき補正される（ステップＳ１０５）。より具体的には、（式１）および最適化されたパラメータセット（ａ，ｂ，ｃ，ｄ，ｅ，ｆ，ｇ）により表される背景モデルが有する輝度プロファイルを有するモデル画像が仮想的に作成される。モデル画像を構成する各画素は、その座標位置（ｘ，ｙ）を最適化された（式１）に代入することで得られるＦ（ｘ，ｙ）の値を輝度値として有する。該モデル画像を基準画像と見なしてシェーディング補正処理を実行することで、処理対象画像における明暗ムラが解消された補正画像が作成される。

この場合の補正処理は、公知のシェーディング補正処理と同様の原理に基づいて行うことが可能であり、例えば次のようにすることができる。すなわち、処理対象画素中の任意の注目画素の輝度値をＬ１、最適化された背景モデルで表されるモデル画像中で、注目画素の座標位置に対応する座標位置にある画素の輝度値をＬ２としたとき、補正画像において当該座標位置を占める画素の輝度値Ｌ３を次式：
Ｌ３＝Ａ・（Ｌ１／Ｌ２）＋Ｂ … （式６）
により求めることで、補正後の画素の輝度値が決定される。

（式６）における符号Ａ，Ｂは予め適宜に定められた実定数であり、このうち定数Ａは０でないものとする。このようにして輝度値が求められた各画素をそれぞれ所定の座標位置に配置することで、原画像における明暗ムラが解消された補正画像が得られる。

なお、この処理において補正の対象となる処理対象画像は、ステップＳ１０３においてウェル内の不要オブジェクトが除去されたものである。気泡や異物など明らかな不要物のみが除去された画像はそのまま補正処理の対象に供することができる。一方、単に背景モデルの作成を目的とするだけであれば、本来の撮像対象物であるウェル内の細胞等に対応するオブジェクトの少なくとも一部も除去された画像や、画像の濃度変化を平滑化するフィルタリング処理が施された画像などを処理対象画像として用いてモデル化処理を実行することが可能である。

このような場合には、撮像対象物に関する画像情報の一部が失われているおそれがあるため、モデル化処理に使用された画像を補正処理の対象とすることは適当でない。したがって、ステップＳ１０１またはＳ１０２で得られた原画像をそのまま補正処理の対象とするか、あるいは当該画像から明らかな不要オブジェクトのみが除去された画像を処理対象画像として別途準備し補正処理に供することが望ましい。

図８は補正画像の例を示す図である。より具体的には、図８（ａ）は、図６（ａ）に示す処理対象画像をこの実施形態のシェーディング補正処理により補正して得られた補正画像Ｉｃの例を示す図であり、図８（ｂ）は当該画像Ｉｃを構成する画素の輝度値の分布を表すヒストグラムである。図６（ａ）に示す処理対象画像Ｉｐと比較すると、図８（ａ）に示す補正画像ＩｃではウェルＷ内の全体がより明るくなっており、またウェル中央部と周縁部との明るさの差も小さくなっている。そして、図８（ｂ）に示すように、各画素の輝度値のばらつきが小さくなっており、明暗ムラが解消されて画像の明るさが均一化されていることがわかる。

撮像対象物としての細胞等を含む画像に対して上記処理が実行されることにより、背景が略均一化されたウェル領域内に細胞等が分布した状態の画像が、補正画像として得られることになる。これにより、例えば液面のメニスカスに起因してウェルＷ内の照明条件が不均一となった原画像からでも、明るさのムラを補正して撮像対象物が明瞭となった画像を作成することが可能となる。

以上説明したように、上記実施形態においては、図４のステップＳ１０１、Ｓ１０３、Ｓ１０４、Ｓ１０５がそれぞれ本発明の「画像取得工程」、「オブジェクト除去工程」、「近似工程」、「補正工程」に相当している。

また、上記実施形態では、ウェルプレートＷＰが本発明の「試料容器」に相当しており、各ウェルＷが本発明の「凹部」に相当している。また、撮像された画像におけるＸ座標軸およびＹ座標軸がそれぞれ、本発明の「第１の座標軸」および「第２の座標軸」に相当している。また、図７（ａ）の二値化画像Ｉｂにおける高輝度領域Ｒｈが、本発明の「主領域」に相当している。

なお、本発明は上記した実施形態に限定されるものではなく、その趣旨を逸脱しない限りにおいて上述したもの以外に種々の変更を行うことが可能である。例えば、上記実施形態は原画像から取得された情報に基づき作成した背景モデルを用いて原画像を補正する機能を有しているが、背景モデルを作成する機能と、それに基づき補正を実行する機能とが分離されてもよい。例えばシェーディング補正機能を担う外部装置等に対して補正に必要なデータを供給するという機能のみを有する装置にも、本発明を適用することが可能である。

この場合、背景モデルを表す例えば（式１）に関するデータ、または該背景モデルで表される仮想的なモデル画像に対応する画像データが本発明の「シェーディング基準データ」に相当する。そして、当該装置が本発明にかかる「シェーディング基準データ作成方法」を実施する主体となる。このような態様では、ステップＳ１０４は本発明の「設定工程」に相当するものとなる。

また、上記実施形態では１回の撮像により１つのウェルＷ全体が撮像される。一方、より大径のウェル（例えば、１つのウェルプレートにおける配設数が６〜４８程度であるウェル）やディッシュなどを撮像する場合には、撮像されるべき範囲が撮像部１３の視野より大きくなる場合がある。この場合には、被撮像範囲を複数の部分画像に分割して撮像し、それらの部分画像を画像処理により結合して被撮像範囲全体をカバーする画像を作成することになる。このような場合でも、個々の部分画像に対して本発明にかかる補正処理を適用することが可能である。補正された部分画像を結合することで、画像全体において明暗ムラの少ない画像を得ることができる。

また、上記実施形態は原画像を撮像する撮像機能を備えた装置であるが、本発明の「画像処理方法」および「シェーディング基準データ作成方法」は、自身が撮像機能を備えず、外部装置等から与えられる原画像データに基づいて処理を行う装置によって実現される態様とすることも可能である。

また、上記実施形態のモデル化処理では、背景モデルを二次元正規分布の確率密度関数によって表すこととしているが、図５に示すモデル化処理自体は、任意の数式で表される背景モデルに対して同じように適用可能である。例えば、二次元正規分布を含む複数種の背景モデルを予め用意しておき、ステップＳ２０１においていずれかの背景モデルが選択される構成とすることが可能である。

また、上記実施形態の補正処理（図４のステップＳ１０５）では、パラメータが最適化された（式１）で表される背景モデルからモデル画像を作成するとして説明しているが、処理対象画像の各画素に対しその輝度値を補正するために必要な情報が提供されれば足り、仮想的なモデル画像の画像データを作成する必要は必ずしもない。例えば、注目画素ごとにその座標値を（式１）に代入して、モデル画像における当該画素の輝度値を都度算出するようにしてもよい。また、最適化された（式１）に座標値を代入して得られる値が（式６）のＬ２に相当するから、例えば（式６）のＬ２に（式１）の右辺を代入して得られた計算式を用いて、注目画素の輝度値と座標値とから補正後の画素値が直接算出されるようにしてもよい。

また、上記実施形態では（式６）を用いて各画素の補正後の輝度値を算出しているが、この方法は一例にすぎず、背景モデルが特定された後の補正処理の具体的内容は上記に限定されるものではなく任意である。

本発明は、撮像により得られた画像、特に、撮像対象物が略一様な背景と共に撮像された画像における明暗ムラを解消することが必要とされる技術分野に好適なものである。

１ 撮像装置
１２ 照明部
１３ 撮像部
１４ 制御部
Ｒｈ 高輝度領域（主領域）
Ｓ１０１ 画像取得工程
Ｓ１０３ オブジェクト除去工程
Ｓ１０４ 近似工程、設定工程
Ｓ１０５ 補正工程
Ｗ ウェル（凹部）
ＷＰ ウェルプレート（試料容器）

Claims (11)

1. 撮像対象物が略一様な背景と共に撮像された原画像を取得する画像取得工程と、
   前記原画像の二次元輝度プロファイルを二次元正規分布の確率密度関数に近似した近似式を特定する近似工程と、
   前記原画像を構成する各画素の輝度値を、前記近似式により表される二次元輝度プロファイルの輝度値に基づいて補正する補正工程と
   を備える画像処理方法。
2. 前記近似工程では、前記原画像の第１の座標軸における座標をｘ、前記第１の座標軸に直交する第２の座標軸における座標をｙとしたときの前記近似式を次式：
   Ｆ（ｘ，ｙ）
   ＝ａ・ｅｘｐ｛−ｂ（ｘ−ｅ）²＋ｃ（ｙ−ｆ）²＋ｄ（ｘ−ｅ）（ｙ−ｆ）｝＋ｇ
   として、未知パラメータａ、ｂ、ｃ、ｄ、ｅ、ｆ、ｇを求めることで前記近似式を特定する請求項１に記載の画像処理方法。
3. 前記近似工程では、Ｌｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔ法により前記近似式を特定する請求項２に記載の画像処理方法。
4. 前記近似工程では、前記原画像の輝度値ヒストグラムにおける最大値、最小値およびモードを求め、
   前記最大値と前記モードとの差の半分の値を閾値として前記原画像を二値化したときに前記閾値以上の輝度値を有する画素が連続してなる領域のうち面積最大のものを主領域としたとき、
   前記未知パラメータそれぞれの初期値を、
   ａ：前記最大値と前記モードとの差、
   ｂ：前記主領域の前記第１の座標軸方向の長さ、
   ｃ：前記主領域の前記第２の座標軸方向の長さ、
   ｄ：０
   ｅ：前記主領域の重心の前記第１の座標軸における座標値
   ｆ：前記主領域の重心の前記第２の座標軸における座標値
   ｇ：前記モード
   とする請求項３に記載の画像処理方法。
5. 前記補正工程では、前記原画像の各画素の輝度値をＬ１、前記近似式により表される二次元輝度プロファイルにおいて対応する画素の輝度値をＬ２としたとき、補正後の当該画素の輝度値Ｌ３を、予め定められた実定数Ａ（Ａ≠０）、Ｂを用いて次式、
   Ｌ３＝Ａ・（Ｌ１／Ｌ２）＋Ｂ
   により求める請求項１ないし４のいずれかに記載の画像処理方法。
6. 前記近似工程に先立って前記原画像中に含まれるオブジェクトの少なくとも１つを前記原画像から除去するオブジェクト除去工程を備える請求項１ないし５のいずれかに記載の画像処理方法。
7. 前記原画像は、液体中の生試料を前記撮像対象物とした画像である請求項１ないし６のいずれかに記載の画像処理方法。
8. 前記原画像は、試料容器に設けられ略円形の水平断面形状を有する凹部に担持された前記液体を前記凹部の全体とともに撮像した画像である請求項７に記載の画像処理方法。
9. 前記原画像は、照明光が前記液体の液面を介して前記撮像対象物に入射した状態で撮像された画像である請求項７または８に記載の画像処理方法。
10. 撮像対象物が略一様な背景と共に撮像された原画像に対しシェーディング補正を行う際の基準となるシェーディング基準データの作成方法において、
    前記原画像を取得する画像取得工程と、
    前記原画像の二次元輝度プロファイルを二次元正規分布の確率密度関数に近似した近似式を特定する近似工程と、
    前記近似式により表される二次元輝度プロファイルを有する仮想的な画像に対応するデータを前記シェーディング基準データとして設定する設定工程と
    を備えるシェーディング基準データ作成方法。
11. 前記近似工程に先立って前記原画像中に含まれるオブジェクトの少なくとも１つを前記画像から除去するオブジェクト除去工程を備える請求項１０に記載のシェーディング基準データ作成方法。