JP2005149500A - 多焦点画像スタック内の歪みの補正方法 - Google Patents

Abstract

【課題】多焦点画像スタック内の歪みを補正する方法を提供する。
【解決手段】多焦点画像スタック内の歪みを補正する方法であって、光学系によって複数のＮのカラー画像を取得するステップであって、各カラー画像は複数の画素を含む、ステップと、Ｎの画像のうちの２個の連続するカラー画像のそれぞれをモノクロ画像に変換するステップと、各モノクロ画像を好都合なＭ×Ｍ画像にスケーリングするステップと、２個の連続する画像のジオメトリック変換およびラジオメトリック変換を計算するステップであって、最小２乗マッチングによって決定される係数の組によってジオメトリック変換およびラジオメトリック変換を決定する、ステップと、２個の連続する画像のそれぞれからの決定済みの係数の組を用いて各画像を順次補正することで、決定されたパラメータの組を連続する各画像に適用することによって多焦点画像を構築するステップと、を含む方法。
【選択図】図７

Description

本発明は、多焦点画像スタック内の歪みを補正する方法に関する。

微細構造の表面を３次元的に再構築するための方法がいくつか知られている。ＣＬＳＭ（共焦点レーザ走査顕微鏡）では、試料を１つの平面内で逐一走査することによって、焦点深度のほとんどない像が得られる。異なる平面での複数の像を用いて、適切な画像処理を用いることで、試料を３次元的に表現できる。ＣＬＳＭでは、データの生成は、非常に費用のかかる方法でないとできず、高品質の光学部品が前提条件である。再構築すべき表面上に細い光線を送ることによって３Ｄ再構築を実現する他の技術においても、高いレベルの技術的に複雑なものが必要である。

非特許文献１の中で言及されている。

特許文献１では、試料の表面プロファイルを光学的に測定するための方法および装置が開示されている。試料の一連のｎ画像が、座標系（ｘ、ｙ、ｚ）のｚ方向の異なる平面内にある画像取得装置によって取得される。結果として生じる画像スタックのｎ画像すべてについて、その画像内容を、ｚ方向の各座標点（ｘ、ｙ）において互いに比較して、そこから所定の基準に従って平面を決定し、その平面の番号（Ｎ）をその座標点（ｘ、ｙ）に割り当てて、マスク画像に記憶する。マスク画像には、試料表面の３Ｄデータがすべて含まれている。２Ｄ画像処理方法を用いて、処理を行なうことができる。マスク画像から、３Ｄ情報を迅速かつ容易に取り出すことができる。表面プロファイルを、３次元的に再構築して表示することができる。開示されたシステムには、多焦点画像スタック内の歪みを補正することは、まったく含まれていない。

特許文献２では、画像処理を用いて１つまたは複数の焦点画像における空間的な関係を変えることで、画像間のジオメトリックな対応関係を保証している。画像処理は、２つの基本プロセス、マッピングおよびリサンプリングからなる。マッピングでは、入力画像の座標点と処理された（出力）画像の座標点との間の空間的な関係を確立する。リサンプリングでは、入力画像からの画素から補間された画像値を、マッピングに従って出力画像の対応する位置に書き込む。隣接する画像をマッチングさせることで補正を決定する提案される方法と比較して、先行技術ではジオメトリの正確な知識に基づいている。加えて、層ごとのラジオメトリック変動に対する補正については、開示もまたは言及もしていない。

特許文献３では、画像処理能力を高めるために複数の画像を用いる顕微鏡システムが開示されている。顕微鏡システムは、ステージを含む顕微鏡、少なくとも１つの拡大レンズ、レンズ・コントローラ、顕微鏡のステージ上の物体についての複数の画像を捕らえるための顕微鏡に結合されたビデオ・キャプチャ・デバイス、複数の画像をビデオ・キャプチャ・デバイスから受け取る処理サブシステムを含む。処理サブシステムは、結果として少なくとも１つの画像を、複数の画像の関数として生成する。処理サブシステムによって生成される、結果としての少なくとも１つの画像には、モザイク、サブモザイク、または一連のモザイクが含まれ得る。システムは、多焦点画像を撮ることはできない。

ＤＥ１０１４９３５７ 米国特許第６，１４８，１２０号明細書 米国特許第６，３１３，４５２号明細書 ＤＥ１０１４９３５１ チェンＪ．（Chen J.）、博士論文「３Ｄ測定に対する複数カメラの使用および幾何学的制約（The Use of Multiple Cameras and Geometric Constraints for 3-D Measurement）」（シティ大学（City University）、英国、１９９５年、２月）

本発明の目的は、多焦点画像スタック内の歪みを補正する方法を提供すること、および歪みのないサンプル画像を提供することである。歪みは、物体とレンズとの距離が変化するときに、実体顕微鏡または他の光学系のノン・テレセントリック・オプティックスから生じる。

上記目的は、請求項１に開示されるステップを含む歪みの補正方法によって達成される。

有利であるのは、本発明の方法が、第１のステップでは、光学系によって複数のＮのカラー画像を取得することである。各カラー画像は、複数の画素を含み、その前の画像または次の画像とは、焦点位置が異なっている。第２のステップでは、Ｎの画像のうちの２個の連続するカラー画像のそれぞれをモノクロ画像に変換することを行なう。第３のステップでは、各モノクロ画像を好都合なＭ×Ｍ画像にスケーリングすることを行なう。第４のステップおよび第５のステップは並行して行なう。第４のステップでは、２個の連続する画像のジオメトリック変換を計算し、第５のステップでは、２個の連続する画像のラジオメトリック変換を計算する。ジオメトリック変換およびラジオメトリック変換を、最小２乗マッチングによって決定される係数の組によって決定する。これは、第６のステップで行なう。第７のステップでは、２個の連続する画像のそれぞれからの決定済みの係数の組を用いて各画像を順次補正することで、決定されたパラメータの組を連続する各画像に適用することによって多焦点画像を構築する。

本発明を、図面に示した典型的な実施形態を参照して、以下に詳細に説明する。

図１は、多焦点画像スタックを取得するための光学系２を示す概略図である。ここで扱う実施形態は、顕微鏡１を示し、これに、ディスプレイ６および入力デバイス８付きのコンピュータ４と、さらに制御ユニット１０とが設けられている。制御ユニット１０は、種々の顕微鏡機能を制御およびモニタするためのものである。制御ユニットには、メモリ９とマイクロプロセッサ１１とが含まれている。メモリ９とマイクロプロセッサ１１とは、多焦点画像スタック内の歪みを補正するための方法によって用いられる。多焦点画像スタック内の歪みを補正する当業者にとっては、顕微鏡は、倒立顕微鏡の可能性もあるしまたは実体顕微鏡の可能性もあることが明らかである。図１に示す光学系２を示すグラフによって、本発明の範囲を限定してはならない。顕微鏡１は、スタンド１２を有し、これに、少なくとも１つの接眼レンズ１４、少なくとも１つの対物レンズ１６、および顕微鏡ステージ１８が取り付けられている。図１に示す実施形態は、顕微鏡ステージ１８を有し、これは直交座標系の３方向に移動させることができる。多焦点画像スタックを得るためには、ステージ１８上のサンプルと対物レンズ１６との間で、相対的な動作を起こさなければならない。相対的な動作は、光学系２の光軸の方向である。ここで示す実施形態では、顕微鏡１はターレット１５を有し、これに少なくとも１つの対物レンズ１６が取り付けられている。少なくとも１つの対物レンズ１６は動作位置にあり、そのため光軸１３（破線で示す）が規定されている。加えて、フォーカシング・ノブ２０が、スタンド１２の両側に設けられている。フォーカシング・ノブ２０を用いて、対物レンズ１６に対する顕微鏡ステージ１８の高さを調節する。顕微鏡１の顕微鏡ステージ１８は、第１のモータ２１によってＸ軸Ｘに沿って動くことができ、第２のモータ２２によってＹ軸Ｙに沿って動くことができ、第３のモータ２３によってＺ軸Ｚに沿って動くことができる。第１のモータ２１、第２のモータ２２、および第３のモータ２３は、制御ユニット１０によって作動される。カメラ２５が、顕微鏡１に取り付けられている。カメラによって、サンプルの画像が記録される。サンプルの画像は、動作位置にある対物レンズ１６によって取られる。電気ケーブル２６によって、カメラと制御ユニット１０とが接続されている。制御ユニット１０は、少なくとも１本の第２の電気ケーブル２７によって、顕微鏡１と接続されている。カメラ２５は、ビデオ・カメラまたはＣＣＤカメラとして構成されていても良い。通常は、図１に示すような複合顕微鏡のビーム経路は、テレセントリックである。すなわち焦点面での倍率と焦点がずれた面での倍率とが同じ値に保たれる。それでも、顕微鏡はテレセントリックの前提条件を満たすため、この場合にも本発明は適用される。

図２に、実体顕微鏡３０の光軸１３とは平行でない軸３３に沿って平行移動した第１の位置３１および第２の位置３２にある実体顕微鏡３０を示す部分図の配置を示す。明瞭にするために、実体顕微鏡３０のいくつかの部分が図示されていない。その部分はたとえば、スタンド１２、制御ユニット１０、カメラ２５、ステージ１８などである。図２では、破線が、対物レンズ１６の第１の位置３１を示している。実線が、対物レンズの第２の位置３２を示している。第１の位置３１から第２の位置３２へ平行移動する結果、実体顕微鏡の光軸１３が、第１の位置１３_１から第２の位置１３_２へ移動して、焦点の合った領域３５_１および３５_２のそれぞれにおいてサンプルから画像が取得される。結果として得られる第１の画像３６_１は、結果として得られる第２の画像３６_２と比べて位置がずれている。位置のずれは、軸３３に沿っての平行移動とともに機械的および光学的効果から生じる他の歪みに起因する。本発明の目的は、一連の画像における連続する画像対の間の歪みを、連続する画像対を適切な配列へ戻す変換のパラメータを見つけることによって、補正することである。

実体顕微鏡では、ノン・テレセントリック・オプティクスが用いられる。このデザインの結果のうち、多焦点画像または画像スタックの忠実度に影響するものがいくつかある。また複合顕微鏡によって生じる多焦点画像または画像スタックの忠実度に影響する問題がいくつか存在する。

連続する画像の横方向の位置ずれを補正するためには、横方向の微調整が必要である。顕微鏡または実体顕微鏡が機械的に不完全であることによって、画像スタック内の単一の画像の横方向の位置ずれが生じること、または画像スタックを通して系統的に進行する位置ずれが導入されることがあり得る。したがって、自動補正または最小２乗マッチングによって横方向の調整をチェックすること、および最終的な画像または画像スタックを生成する前に画像を並べることが必要となる場合があり得る。

集光角度によって、非対称な焦点ずれのにじみが生じる。自動化された多焦点によって、観察ビーム経路の軸と焦点駆動の移動軸とが等しいすべての顕微鏡について、被写界深度が増した優れた画像が生成される。このことは、複合顕微鏡には当てはまるが、実体顕微鏡の場合には、ＡＸ顕微鏡キャリアが取り付けられたものでなければ当てはまらない。したがって、非対称な焦点ずれのにじみを減らすためには、実体顕微鏡の集光角度も考慮することが好都合である。この場合に、ユーザは集光角度と観察ビーム経路の方向とを、これらが自動補正または最小２乗マッチングによって自動的に検出されない場合に、入力できなければならない。集光角度および／または自動的に検出された集光角度を、制御ユニット１０へ入力する。

テレセントリックおよびノン・テレセントリックのビーム経路も補正が必要となる場合がある。通常は、複合顕微鏡のビーム経路はテレセントリックである。すなわち、焦点面での倍率と焦点がずれた面での倍率とが同じ値に保たれる。このことは、実体顕微鏡には当てはまらない。したがって各単一画像のサイズに対して、それらを結合して最終画像にする前に、微調整を行なう必要がある。サイズ変更の量は、自動補正または最小２乗マッチングによって自動的に見つけなければならない。

これらの必要とされる補正は、多焦点画像を得るために必要である。前述の問題を打開するために、最小２乗マッチング（ＬＳＭ）に基づく順次的なトラッキングが開発されている。ＬＳＭテンプレート・マッチングは、画像スタックからの連続する順次的な画像の対に基づいている。順次的なＬＳＭテンプレート・マッチングを、多焦点画像スタック内の動的で正確な画像位置に対して適用して成功している。

ＬＳＭテンプレート・ベースの画像マッチングは、この実体顕微鏡画像スタックの場合には、より高い位置精度が得られ、複合顕微鏡および実体顕微鏡内のノン・テレセントリック・オプティクスによって引き起こされる種々の誤差が補正される。開発された方法は、２個の連続する画像の画像背景テクスチャが非常に似ている長い順次的な画像スタックに対して特に適している。

以下、処理中の原理的なステップを説明する。簡単にするために、データの不必要な記憶および回収を回避するために中間結果を利用できることを用いることが可能な種々の具体化にて実現できる２つの別個のループとして図示する。

順次的な画像スタック（Ｎの画像）を読み込む。Ｎの画像は、カメラによってリアル・タイムに読み込むか、または制御ユニット１０のメモリ９からもしくはコンピュータ４から回収する。次に妥当性チェックおよびメモリ設定を行なう。第１のループではＮの画像のうちの２個の連続する画像のそれぞれに対して、たとえば次のようにして、カラー画像をモノクロ画像に変換する。
Ｐ_ＲＧＢ（ｘ、ｙ）→Ｐ_Ｒ（ｘ、ｙ）、Ｐ_Ｇ（ｘ、ｙ）、Ｐ_Ｂ（ｘ、ｙ）
Ｐ_ｍｏｎｏ（ｘ、ｙ）＝（Ｐ_Ｒ（ｘ、ｙ）＋Ｐ_Ｇ（ｘ、ｙ）＋Ｐ_Ｂ（ｘ、ｙ））／３

次のステップでは、モノクロ画像Ｐ_ｍｏｎｏを、スケーリング（段落１．１で後述する）および双一次補間（段落１．２で後述する）によって好都合なサイズのＭ×Ｍ画像に変換する。Ｍに対する典型的な値は１２８であることが考えられる。
Ｐ_{ＳｃａｌｉｎｇＤｏｗｎ}（Ｘ_ｎｅｗ、Ｙ_ｎｅｗ）＝Ａ_{ＳｃａｌｉｎｇＤｏｗｎ} Ｐ_ｍｏｎｏ（ｘ、ｙ）
ここで
Ｓ_ｘ＝Ｍ／Ｘ_{ｒｅｓｏｌｕｔｉｏｎ}、Ｓ_ｙ＝Ｍ／Ｙ_{ｒｅｓｏｌｕｔｉｏｎ}

以下のステップでは、２個の連続する画像のジオメトリック変換Ａ_{ｇｅｏｍｅｔｒｙ}を計算しなければならない（段落１．３で後述する）。ここで、
Ａ_{ｒａｄｉｏｍｅｔｒｙ}＝［Ｐ_６、Ｐ_７］

Ｐ_０．．．Ｐ_７の計算については、次のセクションで説明する。

この最小２乗マッチング・パラメータの組を、最小２乗マッチング・パラメータ・スタックに保存する。

図３は、カラー画像４１からの画像内での全体強度のグレースケール表示を生成するプロセスを示すグラフである。カメラ２５が捕らえたカラー画像４１を、赤画像４２、緑画像４３、および青画像４４の３つの成分に分離する。分離された赤画像４２、緑画像４３、および青画像４４のデータを、線形演算４５によって結合して、画像内の全体強度のグレースケール表示４６を生成する。

ループ１の後、画像に対する妥当な範囲を算出する。これは、最小２乗マッチング・パラメータ・スタック（ｘ_ｌｏｗ、ｙ_ｌｏｗ、ｘ_{ｕｐｐｅｒ}、ｙ_{ｕｐｐｅｒ}）内のｘおよびｙシフト・パラメータによって行なう。

ループ２では、Ｎの画像内の各単一画像を計算する。
Ｐ_{Ｒｒｅｓｕｌｔ}（ｘ、ｙ）＝ＡＰ_Ｒ（ｘ、ｙ）
Ｐ_{Ｇｒｅｓｕｌｔ}（ｘ、ｙ）＝ＡＰ_Ｇ（ｘ、ｙ）
Ｐ_{Ｂｒｅｓｕｌｔ}（ｘ、ｙ）＝ＡＰ_Ｂ（ｘ、ｙ）
Ｐ_{Ｒｒｅｓｕｌｔ}（ｘ、ｙ）、Ｐ_{Ｇｒｅｓｕｌｔ}（ｘ、ｙ）、Ｐ_{Ｂｒｅｓｕｌｔ}（ｘ、ｙ）→Ｐ_ＲＧＢ（ｘ、ｙ）

計算のプロセスについては、段落１．４で説明する。
Ａ＝Ａ_{ＳｃａｌｉｎｇＤｏｗｎ}Ａ_{ｇｅｏｍｅｔｒｙ}Ａ_{ＳｃａｌｉｎｇＵｐ}

したがって、

ラジオメトリック変換：
Ａ_{ｒａｄｉｏｍｅｔｒｙ}＝［Ｐ_６、Ｐ_７］

したがって、
Ｐ_ＲＧＢ（ｘ、ｙ）＝Ｐ_６＋Ｐ_７ ^＊Ｐ_ＲＧＢ（ｘ、ｙ）

四辺形（ｘ_ｌｏｗ、ｙ_ｌｏｗ、ｘ_{ｕｐｐｅｒ}、ｙ_{ｕｐｐｅｒ}）に基づく画像を切り取った後に、妥当な正確な画像が得られる。（特許文献４）に開示されるように、多焦点アルゴリズムを用いて最終的な多焦点画像を得る。

提案したプロセスの根本的な理論について以下に説明する。スケーリング、双一次補間、および最小２乗マッチングについて詳細に説明する。

１．１スケーリング
スケーリングは特に、大きなサイズの画像データの場合に必要である。通常は、画像サイズは２０４８×１５３６×２４ビットであり、各画像に対して約９Ｍバイトである。このような大きな画像サイズにＬＳＭテンプレート・マッチを適用することは、計算をする上で非現実的である。したがってスケーリング・ダウンを適用して、大きな画像サイズをより小さい画像サイズに小さくする。このスケーリングの他の結果は、画像内の高周波ノイズが減ることである。

図４は、グレースケール画像４６の画素４６_ｘ、_ｙを、計算効率を目的とした一般的により小さい四角形領域上にマッピングするデータ低減画像５０を表わす図である。マッピングとして、グレースケール画像４６をデータ低減画像５０に変換することを、本質的に、グレースケール画像４６の画素を、計算効率および画像内のノイズ低減を目的とした低減画像５０の一般的により小さい四角形領域上にマッピングするリスケーリング行列５１を見つけることによって行なうものを、計算する。双一次補間によって、グレースケール画像４６内の４つの画素からの寄与を結合して、データ低減画像５０内の各画素を決定する。

スケーリング・ダウンには、２つの部分が含まれる。カラー画像をモノクロのグレースケール画像４６に変換する部分、およびスケーリング・ダウンして均一なモノクロ画像サイズ、通常は１２８×１２８にする部分である（しかし十分な計算資源が与えられていれば、他のサイズも可能である）。この例では、１２８×１２８にスケーリングされた画像を用いる。

カラー画像をモノクロ画像に変換する。１つの方法では、単純に、ＲＧＢ画像のＲ、Ｇ、Ｂ成分の平均を適用して、各画素のグレースケール画像強度値を、次のようにして得る。
Ｐ_ＲＧＢ（ｘ、ｙ）→Ｐ_Ｒ（ｘ、ｙ）、Ｐ_Ｇ（ｘ、ｙ）、Ｐ_Ｂ（ｘ、ｙ）
Ｐ_ｍｏｎｏ（ｘ、ｙ）＝（Ｐ_Ｒ（ｘ、ｙ）＋Ｐ_Ｇ（ｘ、ｙ）＋Ｐ_Ｂ（ｘ、ｙ））／３

カラー画像をモノクロ画像に変換する他の方法も、方法の有効性に影響を与えることなく可能である。たとえば次のこの例では、輝度が維持される。
Ｐ_ｍｏｎｏ（ｘ、ｙ）＝０．３０^＊Ｐ_Ｒ（ｘ、ｙ）＋０．５９^＊Ｐ_Ｇ（ｘ、ｙ）＋０．１１^＊Ｐ_Ｂ（ｘ、ｙ）

それぞれが異なる波長で取得される多くのモノクロ画像からなる画像を、別個の成分の強度を好適に一次結合することによって、同様の方法で処理する。本方法は、任意のこのような結合に対して、および任意のビット深度の画像に対して適用される。

画像を１２８×１２８までスケーリング・ダウンするために、ｘおよびｙ方向の両方に対するスケーリング因子を用いなければならない。ｘおよびｙのスケーリング因子の計算は次の式に基づく。
Ｓ_ｘ＝１２８／Ｘｒｅｓｏｌｕｔｉｏｎ
Ｓ_ｙ＝１２８／Ｙｒｅｓｏｌｕｔｉｏｎ

スケーリング中の精度を維持するために、双一次補間アルゴリズムを適用する。

１．２双一次補間
双一次補間法では、グリッド表面６０を、以下の形式のベクトル方程式によって記述する。
Ｓ（ｖ、ｕ）＝Ｔ_１（１−ｕ）（１−ｖ）＋Ｔ_２（１−ｕ）ｖ＋Ｔ_３ｕ（１−ｖ）＋Ｔ_４ｕｖ
ここで、０≦ｕ、ｖ≦１であり、Ｓ（ｖ、ｕ）はグリッド点であり、Ｐ_１、Ｐ_２、Ｐ_３、Ｐ_４は４つの隣接点であり、これらの隣接点での値はそれぞれＴ_１、Ｔ_２、Ｔ_３、およびＴ_４である。

Ｓ（ｖ、ｕ）は、スケーリングにより結果として生じる画素であり、ｕおよびｖによって表わされる小数の画素サイズである。するべき課題は、４つの隣接点に基づいてｕおよびｖを計算して、補間されたグリッド点の値を計算できるようにすることである。Ｓ（０、０）＝Ｔ_１、Ｓ（１、０）＝Ｔ_２等であることが分かる。またＳ（ｖ、１）は点Ｐ_３とＰ４とを結ぶ線分、Ｓ（０、ｕ）は点Ｐ_１とＰ_３とを結ぶ線分等である。表面が四辺形面であるため、ｕおよびｖについての表面の値を、２次方程式を解くことによって見つけることができる。ベクトル方程式を次のように書き換える。
ｘ＝ｘ_１（１−ｕ）（１−ｖ）＋ｘ_２（１−ｕ）ｖ＋ｘ_３ｕ（１−ｖ）＋ｘ_４ｕｖ
ｙ＝ｙ_１（１−ｕ）（１−ｖ）＋ｙ_２（１−ｕ）ｖ＋ｙ_３ｕ（１−ｖ）＋ｙ_４ｕｖ
ここで、ｘ、ｙはグリッド点の座標であり、ｘ_１．．．ｘ_４、ｙ_１．．．ｙ_４はそれぞれ４つの隣接点の座標である。前記方程式を用いて、結果としての画像内の新しい強度を、データを方程式に適用することで計算することができる。

１．３最小２乗マッチング
図６は、第１の領域６１を第２の領域６４上へマッピングするための一般化された一次変換を示す図であり、第１の領域に、２本の軸についての回転拡大ねじりおよび平行移動を取り入れている。一般化された一次変換６３によって、第１の領域６１を第２の領域６４上にマッピングする。変換には、２本の軸についての回転拡大ねじりおよび平行移動が含まれる。これはジオメトリック変換Ｆ_ｇであり、後で詳細に説明する。同時にラジオメトリック変換Ｇ_ｒを行なって、その結果、スケーリングおよびオフセットによって各画素の輝度に影響を与える。この変換の全体を、８つのパラメータによって記述することができ、これらは本発明によって決定される。それによれば、変換された領域６４の基準領域６５に対する類似性は、最小２乗の反復法を用いて、変換された領域６４と基準領域６５との間のグレースケールの差を最小限にすることによって、最大になる。このようにして、画像面のそれぞれの対の再位置合わせに必要な８つのパラメータを漸次決定することができ、そしてこれらのパラメータを用いて、本来のデータを補正することができる。

実体顕微鏡画像スタックの順次的な画像における画像ごとの歪みの原因となる主な因子は２つある。前述したように、これらはジオメトリック歪みおよびラジオメトリック歪みである。明瞭なマッチングおよび高精度な結果を実現するためには、ジオメトリックおよびラジオメトリック因子の両方とも考慮しなければならない。ジオメトリックおよびラジオメトリック歪みは、画像スタック内に存在するが、未知である。これらの歪みを補正してマッチング精度を高めるために、最小２乗マッチング（ＬＳＭ）法を適用する。ＬＳＭは、あるアレイを、残りのグレー値の差を最小限にする他のアレイへ変換することによって、最適なマッチングを規定するアプローチである。簡単な最小２乗の解によって、２個の画像間の対応する点の対のグレー値の間の残りの差異が最小限になる。その結果、最適なマッチングをもたらす要求される最適な変換が表わされる。この方法では、ジオメトリック差異をアフィン変換によってモデリングできる一方で、ラジオメトリック差異は、付加的および乗法的なパラメータによってモデリングされる。単一点ＬＳＭに対する主な目的は、概略的なマッチング位置を精密化してサブ画素の精度を得ることである。導き出された方法について以下に説明する。

Ｉ_１（ｘ、ｙ）を基準の画像テンプレート、Ｉ_２（ｘ、ｙ）を実際の画像パッチと仮定して、２個の画像間のマッピングを次のように表わすことができる。
Ｉ_２（ｘ_２、ｙ_２）＝Ｇ_ｒ（Ｆ_ｇ（Ｉ_１（ｘ_１、ｙ_２）））

ここでＦ_ｇは、ジオメトリック座標変換、すなわち、次の通りである。
Ｆ_ｇ（ｘ’、ｙ’）＝Ｉ_１（ｘ_１、ｙ_１）

またＧ_ｒは、強度またはラジオメトリック変換、すなわち、次の通りである。
Ｇ_ｒ（ｘ”、ｙ”）＝Ｆ_ｇ（ｘ’、ｙ’）

また実際の画像パッチＩ_２（ｘ、ｙ）は、基準の画像テンプレートＩ_１（ｘ_１、ｙ_１）のジオメトリックおよびラジオメトリック変換である。すなわち、次の通りである。
Ｉ_２（ｘ_２、ｙ_２）＝Ｇ_ｒ（ｘ”、ｙ”）

この基本的なアルゴリズムを、次の２つの一般化されたラジオメトリックおよびジオメトリック方程式において、それぞれ定式化することができる。

ラジオメトリック変換に対しては、次の通りである。
Ｉ_２（ｘ_２、ｙ_２）＝Ｇ_ｒ（ｘ”、ｙ”）

しかし、実際の画像におけるノイズの影響を考慮しなければならず、前述の方程式は次のようになる。
Ｉ_２（ｘ、ｙ）−Ｇ_ｒ（ｘ、ｙ）＝ｎ（ｘ、ｙ）
ここで、
Ｉ_２（ｘ、ｙ）は、実際の画像のグレーレベル関数である。
Ｇ_ｒ（ｘ、ｙ）は、テンプレート画像である。
ｎ（ｘ、ｙ）は、誤差ベクトルである。

したがって、次のようになる。
ここで、
は、ｘ方向の傾きである。また、
は、ｙ方向の傾きである。それぞれ、Ｇ_ｘ、Ｇ_ｙと表わすことができる。

ジオメトリック変換に対しては、次の通りである。
Ｇ_ｒ（ｘ”、ｙ”）＝Ｆ_ｇ（ｘ’、ｙ’）

実際の画像におけるノイズの影響を考慮すると、前述の方程式は次のようになる。
Ｇ_ｒ（ｘ、ｙ）−Ｆ_ｇ（ｘ、ｙ）＝ｖ（ｘ、ｙ）
ここで、
Ｇ_ｒ（ｘ、ｙ）は、実際の画像のジオメトリック・レベル関数である。
Ｆ_ｇ（ｘ、ｙ）は、検索ウィンドウである。
ｖ（ｘ、ｙ）は、誤差ベクトルである。

したがって、
ここで、
は、ｘ方向の傾きである。また、
は、ｙ方向の傾きである。それぞれ、Ｆ_ｘ、Ｆ_ｙと表わすことができる。

ジオメトリック変換は、アフィン変換、たとえば平行移動、回転、およびスケーリングなどによって表わされる。すなわち、以下の通りである。
ｘ’＝ｐ_０＋ｐ_１ｘ_１＋ｐ_２ｙ_１
ｙ’＝ｐ_３＋ｐ_４ｘ_１＋ｐ_５ｙ_１

ラジオメトリック変換は、一次変換、たとえば強度スケーリングおよびシフトなどによって表わされる。すなわち、以下の通りである。
Ｉ_２（ｘ_２、ｙ_２）＝ｐ_６＋ｐ_７ ^＊Ｆ_ｇ（ｘ’、ｙ’）

上記ラジオメトリックおよびジオメトリック変換の両方を、実施の際に組み合わせることで、ラジオメトリック品質およびジオメトリック差異が両方とも考慮されて、前述の２つの方程式を次のように書くことができる。
Ｉ_２（ｘ、ｙ）＝ｐ_６＋ｐ_７ ^＊Ｆ_ｇ（ｐ_０＋ｐ_１ｘ＋ｐ_２ｙ、ｐ_３＋ｐ_４ｘ＋ｐ_５ｙ）

したがって、前記方程式を線形にすることによって、一般化された観察方程式を得ることができる。
Ｉ_２（ｘ_ｉ、ｙ_ｊ）−ｎ（ｘ_ｉ、ｙ_ｊ）＝Ｇ^０（ｘ_ｉ、ｙ_ｊ）＋Ｇ^０ _ｐ６（ｘ_ｉ、ｙ_ｊ）ｄｐ_６＋Ｇ^０ _ｐ７（ｘ_ｉ、ｙ_ｊ）ｄｐ_７
Ｇ_ｒ（ｘ_ｉ、ｙ_ｊ）−ｖ（ｘ_ｉ、ｙ_ｊ）＝Ｆ^０（ｘ_ｉ、ｙ_ｊ）＋Ｆ^０ _ｐ０（ｘ_ｉ、ｙ_ｊ）ｄｐ_０＋Ｆ^０ _ｐ１（ｘ_ｉ、ｙ_ｊ）ｄｐ_１＋Ｆ^０ _ｐ２（ｘ_ｉ、ｙ_ｊ）ｄｐ_２＋Ｆ^０ _ｐ３（ｘ_ｉ、ｙ_ｊ）ｄｐ_３＋Ｆ^０ _ｐ４（ｘ_ｉ、ｙ_ｊ）ｄｐ_４＋Ｆ^０ _ｐ５（ｘ_ｉ、ｙ_ｊ）ｄｐ_５
ここで、
ｐ_０、ｐ_１、ｐ_２、ｐ_３、ｐ_４、ｐ_５は、ジオメトリック変換に対する未知の係数である。
ｐ_６、ｐ_７は、ラジオメトリック補正に対する未知の係数である。
Ｆ^０ _ｐ０（ｘ_ｉ、ｙ_ｊ）＝１
Ｆ^０ _ｐ１（ｘ_ｉ、ｙ_ｊ）＝ｘ_ｉ
Ｆ^０ _ｐ２（ｘ_ｉ、ｙ_ｊ）＝ｙ_ｊ
Ｆ^０ _ｐ３（ｘ_ｉ、ｙ_ｊ）＝１
Ｆ^０ _ｐ４（ｘ_ｉ、ｙ_ｊ）＝ｘ_ｉ
Ｆ^０ _ｐ５（ｘ_ｉ、ｙ_ｊ）＝ｙ_ｊ
Ｇ^０ _ｐ６（ｘ_ｉ、ｙ_ｊ）＝１
Ｇ^０ _ｐ７（ｘ_ｉ、ｙ_ｊ）＝Ｆ_ｇ ^０（ｐ_０＋ｐ_１ｘ_ｉ＋ｐ_２ｙ_ｊ、ｐ_３＋ｐ_４ｘ_ｉ＋ｐ_５ｙ_ｉ）

以下のように仮定する。
Ｐ_Ａ＝［ｐ_０、ｐ_１、ｐ_２、ｐ_３、ｐ_４、ｐ_５、ｐ_６、ｐ_７］^Ｔ
ｘ＝ｐ_０＋ｐ_１Ｘ＋ｐ_２Ｙ
ｙ＝ｐ_３＋ｐ_４Ｘ＋ｐ_５Ｙ
ｉ＝ｐ_６＋ｐ_７Ｉ_１（ｘ、ｙ）

ｘ_１、ｙ_１、およびＩ_１（ｘ_１、ｙ_１）は、基準のテンプレート画像の座標および強度値である。ｘ_２、ｙ_２、およびＩ_２（ｘ_２、ｙ_２）は、実際の画像における理想的な位置および強度に対する座標および強度値である。８つの係数ｐ_０．．．ｐ_７は、最小２乗マッチングによって算出することができる。

Ｐ_Ａ＝（Ａ^ＴＷＡ）^−１Ａ^Ｔｘ
ここで、Ａ^ＴはＡの転置行列であり、３ｎ×８行列である。ｎは、テンプレート画像の全画素数である。

またＷは、重み行列であり、通常は共分散行列として選択される。

ここで、
ｑはグローバル分散である。
ｑ_ｉ（ｉ＝１、２Ｋ、３ｎ）は、ローカル分散である。

この最小２乗マッチングは繰返し法で行なわれ、Ｐ_Ａ ^０の初期値は、次の通りである。

図７に、多焦点画像スタック内の歪みを補正する方法のフローチャートを示す。第１のステップ７１では、光学系によって複数のＮのカラー画像を取得する。各カラー画像は、複数の画素を含み、その前の画像または次の画像とは、焦点位置が異なっている。第２のステップ７２では、Ｎの画像のうちの２個の連続するカラー画像のそれぞれをモノクロ画像に変換することを行なう。第３のステップ７３では、各モノクロ画像を好都合なＭ×Ｍ画像にスケーリングすることを行なう。第４のステップ７４および第５のステップ７５は並行して行なう。第４のステップ７４では、２個の連続する画像のジオメトリック変換を計算し、第５のステップ７５では、２個の連続する画像のラジオメトリック変換を計算する。ジオメトリック変換およびラジオメトリック変換を、最小２乗マッチングによって決定される係数の組によって決定する。これは、第６のステップ７６で行なう。第７のステップ７７では、２個の連続する画像のそれぞれからの決定済みの係数の組を用いて各画像を順次補正することで、決定されたパラメータの組を連続する各画像に適用することによって多焦点画像を構築する。

決定された係数の組によって、画像のエッジを切り取る。これらの画像エッジは、得られた決定済みの係数に基づいて妥当でない画像データを有する。画像のスケーリング（第３のステップ７３）は、グレースケール画像４６内の４つの画素からの寄与を結合して、データ低減画像５０内の各画素を決定する双一次補間である。

多焦点画像スタックを取得するための光学系を示す概略図である。 実体顕微鏡の光軸とは平行でない軸に沿って平行移動した第１の位置および第２の位置にある実体顕微鏡を示す部分図の配置である。 カラー画像からの画像内の全体強度のグレースケール表示を生成するプロセスを示すグラフである。 グレースケール画像の画素を、計算効率を目的とした一般的により小さい四角形領域上にマッピングするデータ低減画像を表わす図である。 双一次補間法を示すグラフである。 第１の領域を第２の領域上へマッピングするための一般化された一次変換を示す図であり、第１の領域に、２本の軸についての回転拡大ねじりおよび平行移動を取り入れている。 本発明のプロセスを示すフローチャートである。

符号の説明

１ 顕微鏡
２ 光学系
４ コンピュータ
６ ディスプレイ
８ 入力デバイス
９ メモリ
１０ 制御ユニット
１１ マイクロプロセッサ
１２ スタンド
１３ 光軸
１４ 接眼レンズ
１５ ターレット
１６ 対物レンズ
１８ 顕微鏡ステージ
２０ フォーカシング・ノブ
２１ 第１のモータ
２２ 第２のモータ
２３ 第３のモータ
２５ カメラ
２６、２７ 電気ケーブル
３０ 実体顕微鏡
３１ 第１の位置
３２ 第２の位置
３３ 軸
３５ 焦点の合った領域
３６ 画像
４１ カラー画像
４２ 赤画像
４３ 緑画像
４４ 青画像
４５ 線形演算
４６ グレースケール画像
５０ データ低減画像
５１ リスケーリング行列
６０ グリッド表面
６１ 第１の領域
６３ 一次変換
６４ 第２の領域
６５ 基準領域
７１ 第１のステップ
７２ 第２のステップ
７３ 第３のステップ
７４ 第４のステップ
７５ 第５のステップ
７６ 第６のステップ
７７ 第７のステップ

Claims (11)

1. 多焦点画像スタックにおける歪みを補正する方法であって、
   −光学系によって複数のＮのカラー画像を取得するステップであって、各カラー画像は複数の画素を含むようなステップと、
   −Ｎの画像のうちの２個の連続するカラー画像のそれぞれをモノクロ画像に変換するステップと、
   −各モノクロ画像を好都合なＭ×Ｍ画像にスケーリングするステップと、
   −上記２個の連続する画像のジオメトリック変換およびラジオメトリック変換を計算するステップであって、最小２乗マッチングによって決定される係数の組によって上記ジオメトリック変換およびラジオメトリック変換を決定するようなステップと、
   −上記２個の連続する画像のそれぞれからの決定済み係数の組を用いて各画像を順次補正し、それによって決定されたパラメータの組を連続する各画像に適用することによって多焦点画像を構築するステップと、
   を含む方法。
2. 光学系による複数のＮのカラー画像の取得を、その前の画像および次の画像と焦点位置が異なるように行なう、請求項１に記載の方法。
3. 決定された係数の組に従って画像のエッジを切り取り、画像は、得られた決定済みの係数に基づく妥当でない画像データを有する、請求項１に記載の方法。
4. 画像のスケーリングは、グレースケール画像内の４つの画素からの寄与を結合してデータ低減画像内の各画素を決定する双一次補間である、請求項１に記載の方法。
5. スケーリングされた画像は、通常１２８×１２８ビットである標準的な画像サイズを有する、請求項４に記載の方法。
6. ジオメトリック変換を、８つのジオメトリック係数に対する最小２乗テンプレート・マッチングによって行ない、ラジオメトリック変換を、２つのラジオメトリック係数に対する最小２乗テンプレート・マッチングによって行なう、請求項１に記載の方法。
7. 複数の順次的なＮの画像をカメラによってリアル・タイムに読み込むか、または制御ユニットのメモリから若しくはコンピュータから回収するとともに、妥当性チェックおよびメモリ設定を行なう、請求項１に記載の方法。
8. ジオメトリック変換は、２本の軸についての回転拡大ねじりおよび平行移動を含み、ラジオメトリック変換は、スケーリングおよびオフセットによって各画素の輝度に影響を与える、請求項１に記載の方法。
9. ジオメトリック変換およびラジオメトリック変換を同時に行なう、請求項８に記載の方法。
10. ジオメトリック変換およびラジオメトリック変換の両方を結合して１つの最小２乗マッチング計算にして、一般化された方式を生成する、請求項１〜９のいずれか一項に記載の方法。
11. ジオメトリック変換およびラジオメトリック変換を、２個の連続する画像の対のそれぞれに順次適用して、順次的な画像マッピングを行なう、請求項１０に記載の方法。