JP7748757B2

JP7748757B2 - 顕微鏡画像情報処理方法、顕微鏡画像情報処理システム、およびコンピュータプログラム

Info

Publication number: JP7748757B2
Application number: JP2024517890A
Authority: JP
Inventors: 宣仁森; 泰之木田
Original assignee: National Institute of Advanced Industrial Science and Technology AIST
Current assignee: National Institute of Advanced Industrial Science and Technology AIST
Priority date: 2022-04-26
Filing date: 2023-03-10
Publication date: 2025-10-03
Anticipated expiration: 2043-03-10
Also published as: WO2023210185A1; JPWO2023210185A1

Description

本発明は、顕微鏡画像情報を処理する方法に関する。

従来、顕微鏡で観察したスライドガラス標本などを高精細かつ大範囲にデジタル画像化する技術として、バーチャルスライドという技術が存在する。バーチャルスライドは画像データであるため、スライドガラス標本そのものを保管しておくよりも扱いが容易である。バーチャルスライドは、例えば、遠隔病理診断や病理サンプルのデジタル保管などに使用されうる。

バーチャルスライドを生成するためにいくつかの方法が存在し、例えば、Off-line stitchingと、Real time stitching（例えば非特許文献１）という方法が存在する。Off-line stitchingは、スライドガラス標本などの画像を広範囲にわたって複数撮影した後、オフラインで複数の撮影画像をつなぎ合わせて１枚の画像を生成する方法である。Real time stitchingは、スライドガラス標本を観察している最中に、複数の撮影画像をつなぎ合わせて１枚の画像を生成する方法である。また、バーチャルスライドを作成する方法として、Slide scannerといわれる、自動的にスライドガラス標本をスキャンする装置も存在する。

Alessandro Gherardi1 and Alessandro Bevilacqua, "Real-time whole slide mosaicing for non-automated microscopes in histopathology analysis", [online], 2013年3月30日, National Libray of Medicine, [2022年4月18日検索], インターネット<URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3678752/>

しかしながら、Off-line stitchingにおいては、スライドガラス標本の撮影操作中に接合中の画像をプレビューして確認することができないため、いざ接合処理を実行してみると、接合に失敗する、一部画像が欠落する、等の理由により再撮影作業が発生しやすいという問題があった。また、Real time stitchingは、専用のソフトウェアと、当該ソフトウェアに対応した専用のカメラが必要であり、一般的なデジタルカメラを用いることができない。よって、撮影時に臨機応変に撮影条件を変更することは困難であった。また、Slide scannerといわれるバーチャルスライドを生成するための自動化された専用の装置は、一般的に非常に高価であり、この装置を導入できる環境は限られる。

本発明は、上記課題に鑑みてなされたものである。

上記課題を解決するために、本発明の一態様は、コンピュータシステムによって実行される顕微鏡画像情報処理方法であって、顕微鏡を用いて観察される試料の一部分の撮影画像を取得して記憶領域に保存する第１のステップと、前記撮影画像が前記記憶領域に保存されたことを検知すると、前記撮影画像の特徴点に関する情報である特徴点情報を算出する第２のステップと、前回の前記撮影画像の前記特徴点情報と、前記撮影画像の前記特徴点情報と、を用いて、前回の前記撮影画像の前記特徴点と前記撮影画像の前記特徴点とのマッチング処理を実行する第３のステップと、前記マッチング処理の結果と、現在までに前記記憶領域に保存された複数の前記撮影画像と、に基づいて接合処理を実行して接合済み画像を生成する第４のステップと、前記接合済み画像を表示出力する第５のステップと、を含み、前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップまでを繰り返す、顕微鏡画像情報処理方法である。

また、本発明の他の態様は、顕微鏡を用いて観察される試料の一部分の撮影画像を取得して記憶領域に保存する第１のステップを実行し、前記撮影画像が前記記憶領域に保存されたことを検知すると、前記撮影画像の特徴点に関する情報である特徴点情報を算出する第２のステップを実行し、前回の前記撮影画像の前記特徴点情報と、前記撮影画像の前記特徴点情報と、を用いて、前回の前記撮影画像の前記特徴点と前記撮影画像の前記特徴点とのマッチング処理を実行する第３のステップを実行し、前記マッチング処理の結果と、現在までに前記記憶領域に保存された複数の前記撮影画像と、に基づいて接合処理を実行して接合済み画像を生成する第４のステップを実行し、前記接合済み画像を表示出力する第５のステップを実行し、前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップまでを繰り返す、顕微鏡画像情報処理システムである。

また、本発明の他の態様は、上記の顕微鏡画像情報処理方法をコンピュータシステムに実行させる、コンピュータプログラムである。

本発明の一実施形態に係る顕微鏡画像情報処理方法における逐次接合処理の概要を説明する図である。本発明の一実施形態に係る顕微鏡画像情報処理方法における全体構成処理の概要を説明する図である。Ｎ番目の画像において検出された複数の特徴点と、直前の接合済み画像において検出された複数の特徴点とのマッチングを行った結果の一例を示す図である。連結無向グラフを用いたグローバル変換行列の再計算処理について概念的に説明する図である。連結無向グラフを用いたグローバル変換行列の再計算処理について概念的に説明する図である。本発明の一実施形態に係る顕微鏡画像情報処理方法における逐次接合処理のフローチャートの一例を示す図である。逐次接合処理により生成された接合済み画像（一部）の一例（歪み補正無し）を示す図である逐次接合処理により生成された接合済み画像（一部）の一例（歪み補正あり）を示す図であるステップＳ１２０（２回目以降の画像ファイルに対する処理）における処理のフローチャートの一例を示す図である。本発明の一実施形態に係る顕微鏡画像情報処理方法における全体構成処理のフローチャートの一例を示す図である。全体構成処理のバンドル調整（歪み補正）を行った接合済み画像（一部）の一例を示す図である。全画像を接合して生成される画像の一例を示す図である。矩形領域が複数のタイルに分割される一例を示す図である。逐次接合処理における、領域重複画像の追加マッチング処理の効果を示す図である。逐次接合処理における、領域重複画像の追加マッチング処理の効果を示す図である。逐次接合処理における、Max Spanning Treeの計算と簡易的なバンドル調整による画像どうしのズレの軽減効果を示す図である。逐次接合処理における、Max Spanning Treeの計算と簡易的なバンドル調整による画像どうしのズレの軽減効果を示す図である。全体構成処理における歪み補正の効果を示す図である。全体構成処理における歪み補正の効果を示す図である。コンピュータ装置３０のハードウェア構成の一例を示す図である。

以下、図面を参照しながら本発明の実施形態について詳しく説明する。
（顕微鏡画像情報処理方法の概要）
図１および図２は、本実施形態に係る顕微鏡画像情報処理方法の概要を説明する図である。本実施形態に係る顕微鏡画像情報処理方法は、既存の一般的な顕微鏡１０、カメラ２０、およびコンピュータ装置３０を含んで構成される顕微鏡画像情報処理システム１によって実行されうる。本実施形態に係る顕微鏡画像情報処理方法は、大きく分けて、逐次接合処理と全体構成処理の２つのフレーズによって構成される。逐次接合処理は、スライドガラス標本の一部の画像が撮影されるたびに最新の接合済み画像を生成し、最終的にはスライドガラス標本全体の接合済み画像を生成する処理である。また、全体構成処理は、スライドガラス標本の撮影完了後に実行される処理であり、逐次接合処理が施された接合済み画像に対して実行される処理である。全体構成処理は、主に接合済み画像の質を向上させる各種処理と処理後の画像の保存を行う。逐次接合処理では主に少ない待ち時間で処理が実行され、全体構成処理では主に時間または計算コストのかかる処理が実行される。

図１は、逐次接合処理の概要を説明する図である。ユーザは、顕微鏡１０によってスライドガラス標本を観察しつつ、カメラ用ソフトウェア４０を用いて、スライドガラス標本の各部分をカメラ２０により順次撮影していく。撮影されたスライドガラス標本の各部分の画像５０は、撮影のたびに、カメラ用ソフトウェア４０によって、コンピュータ装置３０が備えるハードディスク装置（ＨＤＤ）３２に書き込まれる（ＳＴ１）。また、接合ソフトウェア４２によってＨＤＤ３２は監視されており（ＳＴ２）、ＨＤＤ３２に画像ファイル５０の書き込みがあったことが検知されるたびに、書き込まれた画像ファイル５０がＨＤＤ３２から読み出されて（ＳＴ３）、接合処理が実行される（ＳＴ４）。スライドガラス標本の各部分の画像５０が撮影されてＨＤＤ３２に保存されるたびに、ＨＤＤ３２から撮影画像が読み出されて接合処理が実行されるという処理が繰り返されて、撮影のたびに接合済み画像５２が更新されていく。そして、接合済み画像５２が更新されるたびに、更新された接合済み画像５２がディスプレイ等の表示装置によってユーザにプレビュー表示される。本実施形態に係る顕微鏡画像情報処理方法においては、生成途中の接合済み画像５２が随時ユーザにプレビューされるため、ユーザは適切な接合済み画像５２が生成されつつあることを随時確認することが可能である。これにより、従来のOff-line stitchingによるバーチャルスライド生成のように、撮影が終了した後に再撮影作業が発生する可能性を低減させることが可能である。

逐次接合処理は、より具体的には以下のようにして実行されうる。すなわち、前回撮影されてＨＤＤ３２に保存された画像と、新たに撮影されてＨＤＤ３２に保存された撮影画像５０との特徴点のマッチングが実行されて両画像間の変換行列が計算される。さらにその結果を基に、ＨＤＤ３２に保存済みの画像であって移動先領域が重複する他の画像との特徴点のマッチングが実行されて変換行列が計算される。計算対象を保存されている全ての撮影画像とするのではなく、移動先領域が重複する画像のみに限定することで、ユーザの待ち時間を削減することが可能である。また、必要に応じてMax Spanning Treeの計算や簡易的なバンドル調整が行われることにより、画像どうしのズレの蓄積を軽減させることが可能である。

図２は、全体構成処理の概要を説明する図である。例えばユーザによって逐次接合処理が完了した旨の入力がなされると、接合ソフトウェア４２は、画像の質を向上させるための処理を実行し（ＳＴ５）、当該処理実行後の接合済み画像５２をＨＤＤ３２に保存する（ＳＴ６）。

全体構成処理は、より具体的には、例えば以下のように実行されうる。
（１）バンドル調整（全特徴点ペアのズレを最小化するような変換行列と、カメラ２０のレンズの歪み補正の計算）
（２）シームの計算（画像間のベストな切れ目を探索する）
（３）露光補正（画像間において露光時間やヴィネッティング（vignetting）の補正）
（４）ブレンド処理（画像間の切れ目が目立たなくるよう画像を合成する）
（５）コンピュータ装置３０のＨＤＤ３２への全体構成処理済み画像の書き込み
上記処理の一部は全体画像の小領域（タイル）に分割して実施される。この場合、（２）～（５）または（３）～（５）はタイル毎に処理されてよい。これにより一度の処理に使用されるメモリ量や計算時間を削減することが可能である。

なお、本実施形態の説明においては、スライドガラス標本を画像化して生成されるバーチャルスライドを一例として説明するが、これはあくまで一例である。本実施形態に係る顕微鏡画像情報処理方法は、例えば、シャーレで培養している細胞の画像化などにおいても同様に適用可能である。

以下、逐次接合処理および全体構成処理についてさらに詳細に説明する。
（逐次接合処理）
図３を参照し、（Ｎ＋１）回目の逐次接合を行うために（Ｎ：１以上の整数）、Ｎ回目に撮影された撮影画像（以下、説明の便宜上、「Ｌａｓｔ（Ｎ）画像」という）５０´と、新たに撮影されてＨＤＤ３２に保存された（Ｎ＋１）回目の撮影画像（以下、説明の便宜上、「Ｎｅｗ（Ｎ＋１）画像」という）５０とについて、既知の方法を用いて複数の特徴点の検出および特徴量の算出と、検出された特徴点のマッチングが実行される。図３は、Ｌａｓｔ（Ｎ）画像５０´において検出された複数の特徴点と、Ｎｅｗ（Ｎ＋１）画像５０において検出された複数の特徴点とのマッチングを行った結果の一例を示す図である。なお、図３においては説明の便宜上、マッチング結果のうちの一部のみが破線で示されている。次に、当該マッチング結果に基づいて、Ｌａｓｔ（Ｎ）画像５０´とＮｅｗ（Ｎ＋１）画像５０との間で変換行列Ｒ_{Ｎ，Ｎ＋１}が計算される。変換行列Ｒ_{Ｎ，Ｎ＋１}は、Ｌａｓｔ（Ｎ）画像５０´とＮｅｗ（Ｎ＋１）画像５０との間の移動距離および回転量を示す行列（アフィン変換行列）である。すなわち、以下の関係式が成り立つ。

ここで、行列

は、Ｎｅｗ（Ｎ＋１）画像５０の特徴点を表す行列である。また、行列

は、Ｌａｓｔ（Ｎ）画像５０´の特徴点を表す行列である。算出された変換行列Ｒ_{Ｎ，Ｎ＋１}は、配列などの形式でコンピュータ装置３０のＨＤＤ３２等の記憶領域に逐次格納される。また、変換行列Ｒ_{Ｎ，Ｎ＋１}の成分ａ、ｂ、ｃ、ｄ、ｅ、ｆは顕微鏡の性能によって自由度を下げてもよい。例えば、顕微鏡の光学系や稼働部が十分な精度を有する場合、撮影される画像間の関係は平行移動と見なせるため、ａ＝ｅ＝１、ｂ＝ｄ＝０として、ｃ、ｆのみ変数としてよい（自由度２）。また、精度がやや悪い場合やスライドガラス標本の固定が緩い場合など、平行移動だけでなく回転も考慮する必要がある場合は、ａ＝ｃｏｓθ、ｂ＝ｓｉｎθ、ｄ＝－ｓｉｎθ、ｅ＝ｃｏｓθとしてθ、ｃ、ｆを変数としてよい（自由度３）。さらに精度が悪い場合や、撮影途中でレンズ倍率が変更になる場合などは、ａ＝ｔ・ｃｏｓθ、ｂ＝ｓｉｎθ、ｄ＝－ｓｉｎθ、ｅ＝ｔ・ｃｏｓθとしてｔ、θ、ｃ、ｆを変数としてよい（自由度４）。顕微鏡や標本に応じて適切な自由度を設定することにより、接合に要する計算時間やメモリ消費量の節減や、接合の精度向上が期待できる。

次に、Ｎｅｗ（Ｎ＋１）画像５０のグローバル変換行列を計算する。「グローバル変換行列」とは、Ｎｅｗ（Ｎ＋１）画像５０が特定の画像を基準として、当該特定の画像からどれぐらい移動および回転をしたかを表す情報である。ここでは１番目の撮影画像を基準として、Ｎｅｗ（Ｎ＋１）画像５０のグローバル変換行列をＲ_Ｎ＋１とすると、以下の関係が成り立つ。

すなわち、Ｎｅｗ（Ｎ＋１）画像５０のグローバル変換行列Ｒ_Ｎ＋１は、Ｎ回目の逐次接合よりも以前の各接合処理においてそれぞれ算出された、直前の撮影画像と新たな撮影画像との間の各変換行列Ｒ_１，２、Ｒ_２，３、・・・Ｒ_{Ｎ，Ｎ＋１}の積として算出されうる。算出されたＮｅｗ（Ｎ＋１）画像５０のグローバル変換行列Ｒ_Ｎ＋１を用いて、基準となる特定の画像（ここでは１番目の撮影画像）の画像データとＮｅｗ（Ｎ＋１）画像５０の画像データとの接合処理が行われ、接合済み画像５２が生成される。

ところで、上述したような方法によってグローバル変換行列Ｒ_Ｎ＋１を算出することが可能であるが、Ｒ_Ｎ＋１を算出するために複数の行列Ｒ_１，２、Ｒ_２，３、・・・Ｒ_{Ｎ，Ｎ＋１}を掛け合わせていくため、各行列において生じた誤差が蓄積していくことになる。本実施形態においてはこの点を考慮して、さらに以下のような処理を行った。

すなわち、グローバル変換行列Ｒ_Ｎ＋１を算出することにより基準となる画像（１回目の画像）からのＮｅｗ（Ｎ＋１）画像５０の大まかな移動先（以下、「移動先領域」という）が分かる。そして、Ｎｅｗ（Ｎ＋１）画像の移動先領域と、過去に算出された各画像についての移動先領域とを比較することにより、Ｌａｓｔ（Ｎ）画像５０´の他に、Ｎｅｗ（Ｎ＋１）画像５０と移動先領域が重複する過去の接合済み画像が存在するか検索する。ここで、「移動先領域が重複する」とは、例えば、移動先領域間の重複面積が“０”を超える場合に重複すると判定されるようになっていてもよいし、予め定められた値を閾値として移動先領域間の重複面積が当該閾値を超えた場合に重複すると判断されるようになっていてもよい。Ｋ番目（０＜Ｋ＜Ｎ）の接合済み画像５０´がＮｅｗ（Ｎ＋１）画像５０と移動先領域が重複していると判断されたとすると、Ｎｅｗ（Ｎ＋１）画像５０とＫ番目の接合済み画像５０´との間でも特徴点のマッチングが行われ、Ｋ番目の接合済み画像５０´とＮｅｗ（Ｎ＋１）画像５０との間の変換行列Ｒ_{Ｋ，Ｎ＋１}が計算される。前記の一連の処理を「領域重複画像の追加マッチング処理」と呼ぶ。また、これらの情報に基づき、各画像を頂点とした連結無向グラフが計算される。

図４Ａおよび図４Ｂは前記グラフを用いたグローバル変換行列の再計算処理について概念的に説明する図である。グローバル変換行列の再計算前は、図４Ａに示されるように、１番目の画像、２番目の画像、３番目の画像、Ｋ番目の画像、Ｎ番目の画像、Ｎ＋１番目の画像（Ｎｅｗ（Ｎ＋１）画像５０）の順番に各画像を辿るパスが決定されていた。領域重複画像の追加マッチング処理が実行されることにより、例えば、図４Ｂに示されるように、Ｋ番目の画像とＮ＋１番目の画像とがリンクされる。そして、Ｋ番目の画像の特徴点とＮ＋１番目の画像の特徴点とのマッチング処理の結果に基づいて、Ｋ番目の画像とＮ＋１番目の画像との間の変換行列Ｒ_{Ｋ，Ｎ＋１}が計算される。また、例えばMax Spanning Treeを計算することによって全体の画像を辿る最適なパスが決定されうる。Max Spanning Treeは重み付き連結無向グラフにおいて、重み総和が最大となる全域木のことであり、クラスカル法などのアルゴリズムで計算できることが知られている。本発明においては各画像間でマッチングした特徴点の数をグラフの重みとして用いることができるが、別の指標を重みとして用いてもよい。Max Spanning Treeを計算した後、当該Treeの中心を新たな基準とする。この処理により、画像を辿るパスとして図４Ａのように１番目の画像からＮ＋１番目の画像までを順番に辿るパスよりも、全体の中央付近、例えば図４Ｂにおける３番目の画像を基準としたほうが全体のルートが短くなるため（行列を掛け合わせる数が少なくなり誤差が低減される）、３番目の画像が新たな基準の画像として決定されうる。この場合、以下の関係が成り立つ。

上記のようにして計算されるグローバル変換行列Ｒ_Ｎ＋１を用いてＮ＋１番目の撮影画像についての逐次接合処理を実行することにより、より誤差を低減させることができる。
（全体構成処理）
全体構成処理は、スライドガラス標本の撮影が全て完了した後に実行される処理であり、全体構成処理によって接合済み画像の質が向上する。本実施形態における全体構成処理は、より具体的には以下のような各処理が実行される。
（１）バンドル調整（全特徴点ペアのズレを最小化するような変換行列と、カメラ２０のレンズの歪み補正の計算）
（２）シームの計算（画像間のベストな切れ目を探索する）
（３）露光補正（画像間において露光時間やヴィネッティングの補正）
（４）ブレンド処理（画像間の切れ目が目立たなくるよう画像を合成する）
（５）コンピュータ装置３０のＨＤＤ３２への全体構成処理済み画像の書き込み
上記処理の一部は全体画像の小領域（タイル）に分割して実施される。この場合、（２）～（５）または（３）～（５）はタイル毎に処理されてよい。これにより一度の処理に使用されるメモリ量や計算時間を削減することが可能である。
（バンドル調整）
上述したように、逐次接合処理においては特徴点のマッチングが行われる。特徴点が変換行列で示されるように移動および回転した時に、マッチングした特徴点の間でどれぐらいの誤差が出るか（再投影誤差）を全特徴点について計算することによって、全体の誤差が計算されうる。この全体の誤差を最小二乗法を用いて最小化することにより、逐次接合済み画像を全体としてより質の良い画像とすることができる。一例として、レーベンバーグ・マーカート（Levenberg-Marquardt）法を用いることができる。レーベンバーグ・マーカート法は非線形最小二乗問題の解法の一つである。

ここで、本実施形態においては、
ｘ_ｉ：処理時点ｉにおける画像の変換行列の成分
Ｊ_ｉ：処理時点ｉにおけるＪａｃｏｂｉａｎ（右肩のＴは転置を示す）
ｒ_ｉ：処理時点ｉにおける特徴点間の誤差
λ：誤差の大きさに応じて調整されるゼロ以上の実数
Ｉ：単位行列
として式（１）が計算されうる。

上記式（１）により、ある処理時点のパラメータｘに対してパラメータＪとパラメータｒが計算される。また、計算されたパラメータＪとパラメータｒを用いて、次のパラメータｘについてさらにパラメータＪとパラメータｒが計算される。これを繰り返していくことで、マッチングした全特徴点間の再投影誤差を最小化するようなグローバル変換行列を計算し、これによりバンドル調整が実行されうる。
（処理フロー：逐次接合処理）
図５は、逐次接合処理のフローチャートの一例を示す図である。

まず、顕微鏡画像情報処理システム１のユーザから、コンピュータ装置３０上で動作するカメラ用ソフトウェア４０や接合ソフトウェア４２などのアプリケーションを介して、ＨＤＤ３２上の監視対象とするフォルダの指定の入力を受け付ける（ステップＳ１０２）。接合ソフトウェア４２は、一定時間ごとに指定されたフォルダの更新を確認する（ステップＳ１０４）。接合ソフトウェア４が逐次接合処理の終了の旨の入力を受け付けない限り（ステップＳ１０６：Ｎｏ）、逐次接合処理を続行し、ステップＳ１０８に進む。

ステップＳ１０８において、接合ソフトウェア４２がステップＳ１０４にて確認したフォルダについて画像ファイルの更新がなかったと判断した場合には（ステップＳ１０８：Ｎｏ）、ステップＳ１０４に戻る。接合ソフトウェア４２が画像ファイルの更新があったと判断した場合、すなわち新たなスライドガラス標本の画像が撮影されてフォルダに保存された場合には（ステップＳ１０８：Ｙｅｓ）、当該画像ファイルの更新が画像ファイルの初回の保存か、すなわち１回目のスライドガラス標本の撮影画像の保存によるものであるか判断する（ステップＳ１１０）。当該画像ファイルの更新が初回の保存ではない場合には（ステップＳ１１０：Ｎｏ）、２回目以降の保存画像ファイルに対する処理を実行する（ステップＳ１２０）。当該処理について後に詳述する。

画像ファイルの更新が画像ファイルの初回の保存であった場合には（ステップＳ１１０：Ｙｅｓ）、保存された画像ファイルをＨＤＤ３２から読み込み、当該画像ファイルにおける特徴点および特徴量を計算し、当該計算結果を、各画像の特徴点および特徴量の情報を保存するための配列に格納する（ステップＳ１１２）。この時、予め歪みパラメータやヴィネッティングパラメータなどが分かっている場合には、画像ファイルの読み込み時に補正を行ってもよい。カメラ２０の光学系レンズによって撮影画像には歪みが生じうる。一般的に、「歪みのない真の座標ｘ，ｙ，ｚ」から「歪みのある撮影後の座標ｕ，ｖ」の変換は、例えば以下の式で計算されうることが知られている。（例えば、Zhengyou Zhang. A flexible new technique for camera calibration. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 22(11):1330-1334, 2000. https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr98-71.pdf）

ここで、ｋ_ｎ、ｐ_ｎ、ｃ_ｘ、ｃ_ｙ、ｆ_ｘ、ｆ_ｙは歪みパラメータであり、これらの値は歪み補正の実行前に推定しておく。より具体的には、歪み補正の実行前に、同じサイズの正方形が並んだ格子状のパターンを、スライドガラス標本を撮影するカメラ２０と同一のカメラ２０で撮影する。そして、撮影された画像の歪みからパラメータｋ_ｎ、ｐ_ｎ、ｃ_ｘ、ｃ_ｙ、ｆ_ｘ、ｆ_ｙの値が推定されうる。推定されたパラメータの値を用いて、（ｕ，ｖ）→（ｘ，ｙ）という逆変換（ここでは２次元系なのでＺ座標は無視する）をすることで撮影画像の歪みを補正することができる。なお、解析的な逆関数は存在しないため、ニュートン法などの近似解を計算するアルゴリズムで逆変換することができる。また、考慮するのは、一部のパラメータのみでもよい（例えば、パラメータｋ_１の値およびパラメータｐ_１の値のみ考慮し、他のパラメータの値は“０”とみなす、など）。また、他の歪みモデルを用いて歪み補正を行ってもよい。

以上のように撮影画像の歪みを補正し、歪み補正後の各撮影画像を用いて接合済み画像を生成することで、接合済み画像におけるズレを低減させることができる。図６Ａおよび図６Ｂは、本願の発明者らが本実施形態に係る顕微鏡画像情報処理方法を用いて生成した接合済み画像（一部）の一例を示す図である。図６Ａは撮影画像に対して歪み補正を行わずに逐次接合を行った時の接合済み画像を示し、図６Ｂは撮影画像に対して歪み補正を行った後に逐次接合を行った時の接合済み画像を示す。図６Ａの接合済み画像の誤差は約３．３２ｅ＋０５であり、図６Ｂの接合済み画像の誤差は約１．５７ｅ＋０４であった。各撮影画像に歪み補正を施して逐次接合を行うほうが誤差が大きく低減された。なお、誤差（再投影誤差）は、以下の手順で算出された。すなわち、（１）グローバル変換行列を用いて各画像の特徴点の基準座標系（基準となる画像の座標系）における座標を計算する。（２）マッチングした特徴点の間で、基準座標系におけるずれを計算する。理想的にはマッチングした特徴点を基準座標系に移した場合、ぴったりと重なるはずであるが、歪みなどがあるためそうならない。これが誤差となる。

図５に戻り、次に、画像ファイルを予め定められた比率で縮小して、画像データを保存するための配列に格納する（ステップＳ１１４）。縮小画像をＲＡＭ３３あるいはＨＤＤ３２等に保持することで、オリジナルサイズの画像を持つより消費メモリや消費記憶容量を削減することができる。ただし、画像ファイルは縮小せずにフルサイズの画像データが配列に格納されてもよく、これ以降、接合済み画像は縮小画像ではなく、フルサイズの画像として生成され扱われてもよい。

ステップＳ１１４にて縮小された１回目の撮影画像のグローバル変換行列（ここでは単位行列となる）を、グローバル変換行列を保存してくおくための配列に格納する（ステップＳ１１６）。ステップＳ１１４にて保存された縮小画像データとステップＳ１１６にて保存されたグローバル変換行列（単位行列）を用いてプレビュー画像を合成し、プレビュー表示する。なお、１枚目の撮影画像の縮小画像をプレビュー表示してもよい（ステップＳ１１８）。

以上の処理は、ステップＳ１０６において、逐次接合処理が完了した（すなわち、バーチャルスライド生成が完了した）などの理由により、接合ソフトウェア４２が、例えばキーボードやマウス等の入力装置を介してユーザからの逐次接合処理の終了の旨の入力を受け付けるまで（ステップＳ１０６：Ｙｅｓ）、続行する。ユーザからの逐次接合処理の終了の旨の入力を受け付けた場合には、逐次結合処理を終了して全体構成処理に移行する。
（処理フロー：ステップＳ１２０の処理）
図７は、ステップＳ１２０（２回目以降の画像ファイルに対する処理）における処理のフローチャートの一例を示す図である。

まず、ｎ回目（ｎ：２以上の整数）にＨＤＤ３２に保存されたスライドガラス標本の一部の画像である画像ｎを読み込み、特徴点および特徴量を計算する（ステップＳ１２０２）。この時、図５のステップＳ１１２と同様に、予め歪みパラメータやヴィネッティングパラメータなどが分かっている場合には、画像ファイルの読み込み時に補正を行ってもよい。次に、画像ｎと、（ｎ－１）回目に撮影された画像である画像ｎ－１との間で特徴点のマッチングを行い、両画像間の変換行列（以下、「画像ペア行列」という）を計算する（ステップＳ１２０４）。ステップＳ１２０４の処理が成功した場合（ステップＳ１２０６：Ｙｅｓ）、画像ｎの移動先領域を計算し、この計算結果に基づいて、保存済みスライドガラス標本の全画像のうち移動先領域が重複する１または複数の画像（以下「ペア画像」という場合もある。図４Ｂの例における画像Ｋが該当）を選択する（ステップＳ１２０８）。次に、画像ｎとステップＳ１２０８において選択された全てのペア画像との間で特徴点のマッチングを行い、それぞれの画像ペア行列（図４Ｂの例における変換行列Ｒ_{Ｋ，Ｎ＋１}）を計算する（ステップＳ１２１０）。ステップＳ１２１０にて実行される特徴点のマッチングは計算コストの大きい処理であり、保存済み画像の全てについて特徴点のマッチングを行うとユーザの待ち時間が大きくなる。よって、ステップＳ１２０８にて選択された移動先領域が重複する画像のみを対象として特徴点のマッチングを行うことでユーザの待ち時間を削減することが可能である。

次に、画像ｎの特徴点および特徴量を配列に格納する（ステップＳ１２１２）。画像ｎを予め定められた比率で縮小し、縮小された画像の画像データを配列に格納する（ステップＳ１２１４）。ステップＳ１２１０において得られた画像ペア行列を配列に格納する（ステップＳ１２１６）。ステップＳ１２１６において格納された画像ペア行列からグローバル変換行列を計算し、配列に格納する（ステップＳ１２１８）。なお、画像ｎのグローバル行列は、画像ｎおよび画像ｎのペア画像（ステップＳ１２０８にて選択された画像。図４Ｂの例における画像Ｋ）の間の画像ペア行列と、画像ｎのペア画像のグローバル変換行列との積として計算することができる。また、ステップＳ１２１８においては、必要に応じて保存済み画像の全画像のグローバル変換行列の再計算を行うようになっていてもよい。「必要に応じて」とは、例えば、「所定回数ごとに」「ユーザから入力を受け付けた場合」「マッチングした全特徴点間の再投影誤差が一定値に達した場合」等がありうる。また、全画像のグローバル変換行列の再計算を行う場合、特徴マッチング数を重みとしたMax Spanning Treeを計算するようになっていてもよい。Treeの中心の画像（図４Ｂの例における３番目の画像）を新たな基準画像とし（３番目の画像のグローバル変換行列は単位行列となる）、Treeのエッジ順に画像ペア行列の積を計算することで各画像のグローバル行列を計算することができる。また、これらに加えて小回数（例えば２回の繰り返し）のバンドル調整を行い、グローバル変換行列を計算してもよい。バンドル調整は、例えばレーベンバーグ・マーカート法によって、マッチングした全特徴点間の再投影誤差を最小化することによって行うことができる。また、全画像のグローバル変換行列の再計算結果に基づいて、計算済みの画像ペア行列を更新してもよい（例えば、画像Ａおよび画像Ｂのグローバル変換行列をそれぞれＭａ、Ｍｂとすると、画像Ａと画像Ｂとの間の画像ペア行列は行列Ｍａと行列Ｍｂ^－１との積で計算することができる）。これにより次回の再計算時に誤差の少ない状態から再計算を開始することができる。

ステップＳ１２１４にて保存された縮小画像データを用いてプレビュー画像を生成し、ディスプレイ等の表示装置に表示する（ステップＳ１２２０）。ここで、画像ｎ－１までのプレビュー画像は合成済みであるので、この画像ｎ－１のプレビュー画像にサイズが縮小された画像ｎを接合することでプレビュー画像を生成してもよい。また、ステップＳ１２１８において保存済み画像の全画像についてグローバル変換行列の再計算を行った場合は、これまでに生成されたプレビュー画像は破棄し、プレビュー画像を生成してもよい。この後、処理は図５のステップＳ１０４に遷移する。

また、ステップＳ１２０６において、ステップＳ１２０４の処理が成功しなかった場合（ステップＳ１２０６：Ｎｏ）、画像ｎと保存済みの全画像との間で特徴点のマッチングを行い、画像ペア行列を計算する（ステップＳ１２２２）。その結果、少なくとも一組の画像ペア（移動先領域が重複する１または複数の画像のペア）が存在する場合には（ステップＳ１２２４：Ｙｅｓ）、処理はステップＳ１２１２に遷移する。画像ペアが存在しない場合には（ステップＳ１２２４：Ｎｏ）、処理は図５のステップＳ１０４へ遷移する。

なお、ステップＳ１２０４の処理が成功しない場合（ステップＳ１２０６：Ｎｏ）、すなわち、画像ｎと画像ｎ－１との特徴点のマッチングに失敗するケースとは、例えば、スライドガラス標本の撮影を、それまで撮影していた視野とは大きく異なる別の視野から再スタートする場合などに発生しうる（必ずしもスライドガラス標本の全領域を一筆書きでなぞるように撮影するとは限らないため）。また、ステップＳ１２２２の処理は、画像ペア行列が一つでも計算できた時点で打ち切ってもよい。その場合、その後の処理はステップＳ１２０８に移行する。

なお、図５～図７のフローにおいて用いられた、画像ペア行列が保存された配列、縮小された画像が保存された配列、グローバル変換行列が保存された配列、各画像の特徴点および特徴量が保存された配列のデータは、後述する全体構成処理に引き継がれる。また、これらのデータの全てまたは一部をファイルとしてコンピュータ装置３０のＨＤＤ３２等の記憶領域に保存しておくことで、全体構成処理を保留しておき、全体構成処理を後日行う、全体構成処理を別のコンピュータ装置で実行する、等が可能である。
（処理フロー：全体構成処理）
図８は、全体構成処理のフローチャートの一例を示す図である。

逐次接合処理が実行された全画像に対して、バンドル調整が行われる。上述したように、バンドル調整はレーベンバーグ・マーカート法（上記式（１））などによって行われうる（ステップＳ３０２）。ここで、事前にカメラ２０の光学系レンズの歪みパラメータが分かっていない場合、バンドル調整の調整対象パラメータとして、グローバル変換行列に加えて、全画像に共通の歪みパラメータを含め、特徴点座標の歪み補正をした上で計算される再投影誤差を最小化することによって、カメラ２０の光学的な歪みの補正を行うことが可能である。すなわち、この時、上記式（１）のパラメータｘは全画像の変換行列の成分と共通の歪みパラメータとを含む。すなわち、式（１）のパラメータｘ_ｉは、歪みを考慮しない場合、各画像のグローバル変換行列の成分（ａ，ｂ，ｃ，ｄ）が並ぶ（ｘ_ｉ＝（Ｒ_１，Ｒ_２，・・・，Ｒ_ｎ）。ここで、（Ｒ_１＝ａ_１、ｂ_１、ｃ_１，ｄ_１））。歪みを考慮した場合は、ｘ_ｉ＝（Ｒ_１，Ｒ_２，・・・，Ｒ_ｎ，ｋ_１，ｋ_２，ｐ_１，ｐ_２，・・・）のように式（２）の歪みパラメータを追加したかたちになる。これが共通の歪みパラメータになる。このパラメータｘ_ｉを使って誤差ｒを計算すれば、式（１）に基づいて改善されたｘ_{（ｉ＋１）}を計算することができ、これを繰り返すことで、より改善されたグローバル変換行列と歪みパラメータが得られる。初期値としては、歪み無し（ｋ_１＝ｐ_１＝・・・＝０）に仮決めして、式（１）を繰り返し計算することでｋ_１，ｐ_１，・・・が定まっていく。また、各画像に固有の歪みパラメータを想定して計算することも可能だが、共通の歪みパラメータを用いることによって、より計算時間を抑制することが可能であり、スライドガラス標本を観察する場所を変えるだけでカメラ２０のレンズの変更は伴わない場合は充分に実用的である。なお、ステップＳ３０２では特徴点座標のみを用いて歪みパラメータを計算するため、実際に各画像の全画素を対象として歪み補正を行うのは以降の処理である。

次に、縮小画像（ステップＳ１１４およびＳ１２１４において配列に格納された縮小画像）について歪み補正が行われる（ステップＳ３０４）。なお、逐次接合処理に先んじて事前に歪みパラメータが分かっており、図５のステップＳ１１２および図７のステップＳ１２０２において各撮影画像について歪み補正をすでに行っている場合は、本ステップでの補正は省略される。図９は、図６Ａおよび図６Ｂと同じ撮影画像に関して、式（１）のパラメータｘに共通の歪みパラメータを含めることで歪み補正を行ったうえで画像の接合をした場合の接合済み画像の例を示す図である。図６Ａの接合済み画像の誤差は約３．３２ｅ＋０５であり、図６Ｂの接合済み画像の誤差は約１．５７ｅ＋０４であった。図９の接合済み画像の誤差は約１．５５ｅ＋０３であった。図９の接合済み画像の誤差は、図６Ｂに示された、予め撮影画像に対して歪み補正を行って逐次接合を行った接合済み画像よりも、さらに誤差が低減されていることが分かる。

次に、縮小画像を用いて撮影画像間のつなぎ目（シーム）を計算する（ステップＳ３０６）。シームの計算が行われることにより、適切な切れ目においてが画像を接合することができる。シームを計算する方法として、例えば、ボロノイ（Voronoi）図を用いる方法、動的計画法（Dynamic programing）、グラフカット（Graphcut）、などの既存の手法を用いることができる。なお、シームの計算には時間を要するため、本実施形態においては、ユーザの待ち時間が長くとも許容されるように、逐次接合処理ではなく全体構成処理にてシームの計算が実行されることとしている。また、フルサイズの画像を使うとシームの計算に時間がかかるが、本実施形態においては上述したように縮小画像を用いており、これは計算時間の短縮に寄与する。

次に、全撮影画像を接合して生成される画像の大きさ（幅および高さ）が計算される。図１０は、全画像を接合して生成される画像の一例を示す図である。なお、図１０は説明の便宜上、簡潔化して示されている。画像が６枚（画像１～６）あるとすると、画像１～６を接合して生成される画像は、画像１～６全体を含むように構成される画像６０である。本ステップでは、画像６０の大きさ（幅ｗ，高さｈ）を以下のように計算する。すなわち、（１）ステップＳ３０２にて得られたグローバル変換行列と歪み補正パラメータとを用いて、各画像１～６の（基準座標系における）移動先領域を計算し（各画像１～６の大きさ（幅および高さ）から計算することができる）、（２）（１）で得られた全移動先領域を含む矩形領域を計算する。また、計算された領域は複数の小領域（以下、「タイル」「タイル画像」等という）に分割される（ステップＳ３０８）。図１１は、矩形領域６０が複数（ここでは説明の簡便化のため２つ）のタイル６５、６６に分割される一例を示す図である。また、タイル６５、６６とそれぞれ交差する画像をタイル画像セットとする。本例では、タイル６５のタイル画像セットは、画像１、２、３、４、５であり、タイル６６のタイル画像セットは、画像４、５、６である。例えばタイル６５の接合処理をするとタイル６５のサイズより大きい画像が生成されるが、はみ出した部分は切り落とす。また、タイルの大きさや形状は、以下のようにして決定されうる。すなわち、（１）予め定められた大きさの矩形（1024pixel×1024pixelなど）、（２）ユーザに設定ファイルなどで指定させる、（３）コンピュータ装置３０のメモリを考慮して自動決定される、等となっていてよい。（３）については、例えば、「タイル画像セットに含まれる画像の数×フルサイズ画像１枚のサイズ」により１タイルを処理する際に必要なメモリを見積ることができるため、この必要メモリ量がソフトウェアにおいて使用可能なメモリの１０％を超えないようにタイルの大きさを決定する、等となっていてよい。また、タイルの形状は、予め定められた特定の形状（正方形など）としてもよいし、個々の画像１～６の相似形（縦横比３：４の長方形など）としてもよい。以下、ステップＳ３１０～Ｓ３１８の処理では、画像を接合するため、フルサイズの画像が使用される。また、フルサイズの画像が全て一度に処理されるとメモリ使用量と計算時間が大きくなるため、画像全体がタイルに分割されてタイルごとにステップＳ３１０～Ｓ３１８の処理が実行される。

まず、本ステップにおいて処理対象となっているタイルに対して、移動先領域が交差する画像を特定する。図１１の例においては、タイル６５のタイル画像セットは、画像１、２、３、４、５であり、タイル６６のタイル画像セットは、画像４、５、６である。（ステップＳ３１０）。（特定された画像のセットを「タイル画像セット」という。）次に、ステップＳ３１０にて特定されたタイル画像セットの各画像ファイルがＨＤＤ３２から読み出される（ステップＳ３１２）。ステップＳ３１２において読み出されたタイル画像セットについて、歪みパラメータ（予め分かっている歪みパラメータまたは上記式（１）のパラメータｘに共通の歪みパラメータを含ませるようにして推測された歪みパラメータのいずれでもよい）に基づいた歪み補正が行われ、さらに予めヴィネッティングパラメータが分かっている場合には、当該ヴィネッティングパラメータに基づいた画像補正を行う。予めヴィネッティングパラメータが分かっていない場合には、次のステップにて露光補正を行う（ステップＳ３１４）。露光補正の方法は、例えば、M. Uyttendaele, A. Eden and R. Skeliski (2001) "Eliminating ghosting and exposure artifacts in image mosaics," Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001, p. II（https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Uyttendaele01.pdf）に開示された方法を用いることができる。

ステップＳ３０２において算出されたグローバル変換行列に従ってタイル画像セットの全画像を移動して配置する（ステップＳ３１６）。この際、ステップＳ３０６において計算されたシーム情報を用いて全画像を接合する。また、接合の際にブレンド処理を行ってもよい。ブレンド処理は、リニアブレンド、マルチバンドブレンドなど様々な既知の手法が用いられてよい。例えば、Richard Szeliski (2007), "Image Alignment and Stitching: A Tutorial", Foundations and Trends (R) in Computer Graphics and Vision: Vol. 2: No. 1, pp 1-104. （https://www.nowpublishers.com/article/Details/CGV-009）に開示された方法が用いられうる。

ステップＳ３１６において合成された画像（タイル）を一時ファイルとしてコンピュータ装置３０のＨＤＤ３２等の記憶領域に保存する（ステップＳ３１８）。全てのタイルについて、以上のステップＳ３１０～Ｓ３１８の処理が実行される。

各タイルについてステップＳ３１８において保存された一時ファイルをコンピュータ装置３０のＨＤＤ３２等の記憶領域から読み出し、最終ファイルに順次書き込む（ステップＳ３２０）。これにより、全体構成処理済みの最終的なバーチャルスライド画像がコンピュータ装置３０のＨＤＤ３２等の記憶領域に保存される。

図１２Ａおよび図１２Ｂは、逐次接合処理における、領域重複画像の追加マッチング処理の効果を示す図である。図１２Ａのグラフは、横軸が画像の登録順、縦軸がマッチングを試みた画像の数を表す。図１２Ｂのグラフは、横軸が画像の登録順、縦軸がマッチングに要した時間を表す。ここでは逐次接合処理において、三種類の方式を用いた。第一の方式は、新たな画像（Ｎｅｗ（Ｎ+１）画像）を検出した際に、それまでに接合した全ての画像に対してマッチングを試みる方式である。第二の方式は、領域重複画像の追加マッチング処理を行う方式である。第三の方式は、新たな画像（Ｎｅｗ（Ｎ+１）画像）を検出した際に、直前の接合済み画像（Ｌａｓｔ（Ｎ）画像）のみに対してマッチングを試みる方式である。第一の方式では画像の登録順が後になるにつれて登録対象となる画像の数が線形に増加する。このため、マッチングに要する時間も増大しており、ユーザの待ち時間が次第に増加してしまう。第二の方式ではマッチングを試みた画像の数は変動しつつも画像の登録順にかかわらずほぼ一定である。このため、マッチングに要する時間もほぼ定常である。第三の方式ではマッチングを試みる画像の数は常に１であり、従ってマッチングに要する時間も定常である。しかしながら、第三の方式を用いる場合、新たな画像（Ｎｅｗ（Ｎ+１）画像）と、直前の接合済み画像（Ｌａｓｔ（Ｎ）画像）以外の画像との位置関係が分からないため、接合の繰り返しや光学系の歪みに起因するズレを、Max Spanning Treeの計算やバンドル調整により修正することができない。このように、第二の方式（領域重複画像の追加マッチング処理）を用いることで、マッチングに要する時間（ユーザの待ち時間）を定常に保ちつつ、画像どうしのズレを修正することができる。

図１３Ａおよび図１３Ｂは、逐次接合処理における、Max Spanning Treeの計算と簡易的なバンドル調整とによる画像どうしのズレの軽減効果を示す図である。図１３Ａは膵臓がんのＨＥ染色スライド（撮影枚数３２０枚）の場合のグラフを示し、図１３Ｂは肝臓のＨＥ染色スライド(撮影枚数５９枚)の場合のグラフを示す。これらのグラフは、接合を１０回行うごとにMax Spanning Treeの計算と簡易的なバンドル調整とを行った場合と、これらを行わなかった場合とにおける、逐次接合処理が終了した時点での画像どうしのズレの大きさ(ＲＭＳ: root mean square)を示している。いずれのグラフにおいてもMax Spanning Treeの計算と簡易的なバンドル調整とを行った場合において画像どうしのズレが大幅に軽減されていることが分かる。ズレがあまりに大きいとユーザが適切な接合済み画像が生成されているかを確認するのに支障が出るため、必要に応じてMax Spanning Treeの計算や簡易的なバンドル調整によってズレを軽減することは有効である。

図１４Ａおよび図１４Ｂは、全体構成処理における歪み補正の効果を示す図である。これらは、特に、上述したステップＳ３０４の処理に関する。図１４Ａは膵臓がんのＨＥ染色スライド（撮影枚数３２０枚）の場合のグラフを示し、図１４Ｂは肝臓のＨＥ染色スライド(撮影枚数５９枚)の場合のグラフを示す。また、これらの図は、上述した式（１）のパラメータｘに共通の歪みパラメータを含めることで歪み補正を行った場合と、歪み補正を行わなかった場合とにおける、ズレの大きさ(ＲＭＳ)を示している。なお、ここでは、ＲＭＳは特徴点のペア一つあたりのズレをピクセル単位で示した数値である。歪み補正を行った場合、０．５ピクセル以下までズレが低減されていることが分かる。これは目視不能なレベルのズレであり、歪み補正の効果が極めて大きいことを示している。

（顕微鏡画像情報処理システムの構成）
図１および図２に示されるように、本実施形態に係る顕微鏡画像情報処理方法を実行する顕微鏡画像情報処理システム１は、顕微鏡１０と、カメラ２０と、コンピュータ装置３０とを含んで構成されうる。顕微鏡１０、カメラ２０、およびコンピュータ装置３０は、既存の一般的なコンピュータ装置を用いることができる。図１５は、コンピュータ装置３０のハードウェア構成の一例を示す図である。コンピュータ装置３０は、プロセッサ３１、ＨＤＤ３２、ＲＡＭ（Random Access Memory）３３、ＲＯＭ（Read Only Memory）３４、ＣＤ、ＤＶＤ、ＵＳＢメモリ、メモリスティック、ＳＤカード等のリムーバブルメモリ３５、入出力ユーザインタフェース（キーボード、マウス、タッチパネル、スピーカ、マイク、ＬＥＤ（light-emitting diode）等）３６、他のコンピュータ装置と通信可能な有線／無線の通信インタフェース３７、ディスプレイ３８、等の一般的なコンピュータ装置と同様のハードウェア構成を備えうる。コンピュータ装置３０は、例えばＨＤＤ３２に記憶されているコンピュータプログラム（カメラ用ソフトウェア４０、接合ソフトウェア４２、およびその他各種のコンピュータプログラム）および処理対象の各種データをＲＡＭ３３等のメモリに読み出して実行することで、上述した本実施形態に係る顕微鏡画像情報処理方法の各処理を実現しうる。

なお、図１および図２においては、コンピュータ装置３０は１つの装置として図示されているが、２つ以上のコンピュータ装置によって構成されていてもよい。例えば、スライドガラス標本の各部分の画像をＨＤＤ３２に保存する第１のコンピュータ装置と、これとは別体の第２のコンピュータ装置が、第１のコンピュータ装置に保存された複数の撮影画像を用いて逐次接合処理以降の処理を実行するようになっていてもよい。この際、例えば、第１のコンピュータ装置のＨＤＤ３２に一時的に保存されたスライドガラス標本の撮影画像は、自動的に第２のコンピュータ装置へ有線または無線の通信によって送信され、第２のコンピュータ装置はこれをトリガーとして、以降の逐次接合処理および全体構成処理を実行するようになっていてもよい。また、別の方法として、例えば、第１のコンピュータ装置が逐次接合処理を実行し、逐次接合処理が完了した逐次接合済み画像および全体構成処理に必要なデータを第２のコンピュータ装置に有線または無線の通信によって送信するようになっていてもよい。そして、第２のコンピュータ装置は送信された接合済み画像および全体構成処理に必要なデータを用いて全体構成処理を実行するようになっていてもよい。

また、本実施形態においてはスライドガラス標本の各部分の撮影画像や接合済み画像はＨＤＤ３２に保存されることとして説明したが、他の記録媒体に保存されるようになっていてもよい。

本実施形態に係る顕微鏡画像情報処理方法によれば、スライドガラス標本の各部分の画像が撮影されながら逐次接合処理が実行され、接合途中の画像が逐次ユーザにプレビュー表示される。ユーザはプレビュー表示を見て、適切に接合処理が実行されているかをその都度確認しつつ、スライドガラス標本の撮影を進めることができる。これにより、従来のOff-line stitchingのように最終的に接合に失敗して、スライドガラス標本を再撮影する作業が必要となることを回避しうる。

また、本実施形態に係る顕微鏡画像情報処理方法によれば、従来のReal time stitchingのように自動電動ステージや専用のカメラなどは不要であり、既存のイメージングシステム（顕微鏡、カメラ、コンピュータ装置）を使用して実現可能である。これにより、システムの導入コストが低く、日常的に使い慣れた顕微鏡を使用して見慣れた撮影画像でバーチャルスライドを生成可能である。また、バーチャルスライドを作成する従来の方法として、Slide scannerといわれる、自動的にスライドガラス標本をスキャンする装置も存在するが、このような装置は非常に高価であるが、本実施形態に係る顕微鏡画像情報処理方法によれば、このような専用の装置は不要である。

また、蛍光サンプルを撮影するためには露光時間が５００ｍ～１秒程度必要である。また、多重に染色されたサンプルのバーチャルスライドを生成したい場合はフィルタの手動の切り替えが必要となる。例えば従来のReal time stitchingでは、このように撮影条件を臨機応変に変更すること等は困難であった。これに対し、本実施形態に係る顕微鏡画像情報処理方法は、サンプル画像を１枚ずつ撮影していくため蛍光サンプルにも対応可能である。

また、本実施形態に係る顕微鏡画像情報処理方法によれば、プレビュー画像は縮小画像によって生成し、全体構成処理において元のサイズの接合済み画像を生成することができるため、使用するメモリの削減と処理スピードの改善に資する。

ここまで、本発明の一実施形態について説明したが、本発明は上述の実施形態に限定されず、その技術的思想の範囲内において種々異なる形態にて実施されてよいことは言うまでもない。

また、本発明の範囲は、図示され記載された例示的な実施形態に限定されるものではなく、本発明が目的とするものと均等な効果をもたらすすべての実施形態をも含む。さらに、本発明の範囲は、各請求項により画される発明の特徴の組み合わせに限定されるものではなく、すべての開示されたそれぞれの特徴のうち特定の特徴のあらゆる所望する組み合わせによって画されうる。

なお、以下のような構成も本発明の技術的範囲に属する。
（１）
コンピュータシステムによって実行される顕微鏡画像情報処理方法であって、
顕微鏡を用いて観察される試料の一部分の撮影画像を取得して記憶領域に保存する第１のステップと、
前記撮影画像が前記記憶領域に保存されたことを検知すると、前記撮影画像の特徴点に関する情報である特徴点情報を算出する第２のステップと、
前回の前記撮影画像の前記特徴点情報と、前記撮影画像の前記特徴点情報と、を用いて、前回の前記撮影画像の前記特徴点と前記撮影画像の前記特徴点とのマッチング処理を実行する第３のステップと、
前記マッチング処理の結果と、現在までに前記記憶領域に保存された複数の前記撮影画像と、に基づいて接合処理を実行して接合済み画像を生成する第４のステップと、
前記接合済み画像を表示出力する第５のステップと、
を含み、
前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップまでを繰り返す、顕微鏡画像情報処理方法。
（２）
前記第４のステップにおいて生成される接合済み画像は、予め定められた比率で縮小された前記撮影画像を用いて生成される画像である、上記（１）の顕微鏡画像情報処理方法。
（３）
前記第３のステップにおいて、前記記憶領域に保存されている前記撮影画像のうちの特定の画像を基準とした、前記第１のステップにて保存された前記撮影画像のグローバル変換行列が算出される、上記（１）または（２）の顕微鏡画像情報処理方法。
（４）
前記グローバル変換行列はMax Spanning Treeによって算出される、上記（３）の顕微鏡画像情報処理方法。
（５）
前記第３のステップにおける前記マッチング処理が行われた全特徴点間の再投影誤差を最小二乗法を用いて最小化することによって前記グローバル変換行列を調整する、上記（４）の顕微鏡画像情報処理方法。
（６）
前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像を含む領域について、当該領域と、全ての前記撮影画像と、をグローバル変換行列に従って配置するステップを、前記第５のステップの後にさらに含む、上記（１）の顕微鏡画像情報処理方法。
（７）
前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像を含む領域を複数の小領域に分割する第６のステップと、
前記複数の小領域の各々について、
処理対象の小領域と交差する１または複数の前記撮影画像を特定し、
前記処理対象の小領域と、前記特定された１または複数の前記撮影画像と、をグローバル変換行列に従って配置する
ことを実行する第７のステップと、
をさらに含む、上記（１）～（４）の顕微鏡画像情報処理方法。
（８）
前記第６のステップにおいて、前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像に対してバンドル調整を実行してから、接合済み画像を複数の小領域に分割する、上記（７）に記載の顕微鏡画像情報処理方法。
（９）
前記バンドル調整は以下のレーベンバーグ・マーカート法で定義される式を評価することにより実行される、上記（８）に記載の顕微鏡画像情報処理方法。

ただし、
ｘ_ｉ：処理時点ｉにおける画像の変換行列の成分
Ｊ_ｉ：処理時点ｉにおけるＪａｃｏｂｉａｎ（右肩のＴは転置を示す）
ｒ_ｉ：処理時点ｉにおける特徴点間の誤差
λ：誤差の大きさに応じて調整されるゼロ以上の実数
Ｉ：単位行列
（１０）
前記バンドル調整は以下のレーベンバーグ・マーカート法で定義される式を評価することにより実行される、上記（８）に記載の顕微鏡画像情報処理方法。

ただし、
ｘ_ｉ：処理時点ｉにおける画像の変換行列の成分および歪み補正パラメータ
Ｊ_ｉ：処理時点ｉにおけるＪａｃｏｂｉａｎ（右肩のＴは転置を示す）
ｒ_ｉ：処理時点ｉにおける特徴点間の誤差
λ：誤差の大きさに応じて調整されるゼロ以上の実数
Ｉ：単位行列
（１１）
前記試料は蛍光試料である、上記（１）～（７）の顕微鏡画像情報処理方法。
（１２）
前記コンピュータシステムは、第１のコンピュータ装置と、前記第１のコンピュータ装置とは別体である第２のコンピュータ装置と、を含んで構成され、
前記第１のステップは前記第１のコンピュータ装置において実行され、
前記第２のステップ～前記第５のステップは前記第２のコンピュータ装置において実行される、上記（１）～（８）の顕微鏡画像情報処理方法。
（１３）
顕微鏡を用いて観察される試料の一部分の撮影画像を取得して記憶領域に保存する第１のステップを実行し、
前記撮影画像が前記記憶領域に保存されたことを検知すると、前記撮影画像の特徴点に関する情報である特徴点情報を算出する第２のステップを実行し、
前回の前記撮影画像の前記特徴点情報と、前記撮影画像の前記特徴点情報と、を用いて、前回の前記撮影画像の前記特徴点と前記撮影画像の前記特徴点とのマッチング処理を実行する第３のステップを実行し、
前記マッチング処理の結果と、現在までに前記記憶領域に保存された複数の前記撮影画像と、に基づいて接合処理を実行して接合済み画像を生成する第４のステップを実行し、
前記接合済み画像を表示出力する第５のステップを実行し、
前記試料の一部分の撮影が終了するまで、前記第１のステップから前記第５のステップまでを繰り返す、顕微鏡画像情報処理システム。
（１４）
上記（１）～（１２）の顕微鏡画像情報処理方法をコンピュータシステムに実行させる、コンピュータプログラム。

１…顕微鏡画像情報処理システム
１０…顕微鏡
２０…カメラ
３０…コンピュータ装置
３１…プロセッサ
３２…ハードディスク装置（ＨＤＤ）
３３…ＲＡＭ
３４…ＲＯＭ
３５…リムーバブルメモリ
３６…入出力ユーザインタフェース
３７…通信インタフェース
３８…ディスプレイ
４０…カメラ用ソフトウェア
４２…接合ソフトウェア
５０、５０´…スライドガラス標本の撮影画像

Claims

コンピュータシステムによって実行される顕微鏡画像情報処理方法であって、
顕微鏡を用いて観察される試料の一部分の撮影画像を取得して記憶領域に保存する第１のステップと、
前記撮影画像が前記記憶領域に保存されたことを検知したことをトリガーとして、前記撮影画像の特徴点に関する情報である特徴点情報を算出する第２のステップと、
前回までに前記記憶領域に保存された１または複数の前記撮影画像の前記特徴点情報と、前記第２のステップにおいて算出された前記撮影画像の前記特徴点情報と、を用いて、前回までに前記記憶領域に保存された１または複数の前記撮影画像の前記特徴点と前記第１のステップにおいて前記記憶領域に保存された前記撮影画像の前記特徴点とのマッチング処理を実行する第３のステップと、
前記マッチング処理の結果と、前回までに前記記憶領域に保存された１または複数の前記撮影画像と、前記第１のステップにおいて前記記憶領域に保存された前記撮影画像と、を用いて接合処理を実行して接合済み画像を生成する第４のステップと、
前記接合済み画像を表示出力する第５のステップと、
を含み、
前記試料を部分的に撮影するたびに前記第１のステップから前記第５のステップまでを繰り返し、
前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像を含む領域を複数の小領域に分割する第６のステップと、
前記複数の小領域の各々について、
処理対象の小領域と交差する１または複数の前記撮影画像を特定し、
前記処理対象の小領域と、前記特定された１または複数の前記撮影画像と、をグローバル変換行列に従って配置する
ことを実行する第７のステップと、
をさらに含む、顕微鏡画像情報処理方法。
前記第４のステップは、特定の条件下においては、前記マッチング処理の結果と、前回までに前記記憶領域に保存された１または複数の前記撮影画像と、前記第１のステップにおいて前記記憶領域に保存された前記撮影画像と、を用いて接合処理を実行して接合済み画像を生成し、前記特定の条件下ではない場合は、前記マッチング処理の結果と、前回に接合された接合済み画像と、前記第１のステップにおいて前記記憶領域に保存された前記撮影画像と、を用いて接合処理を実行して接合済み画像を生成し、
前記特定の条件は、前記試料の撮影回数があらかじめ定められた回数の倍数であること、ユーザからの入力を受け付けたこと、または前記第３のステップにおいてマッチングされた全特徴点間の再投影誤差が一定値に達したことである、請求項１に記載の顕微鏡画像情報処理方法。
前記第３のステップにおける前記マッチング処理は、
前回までに前記記憶領域に保存された１または複数の前記撮影画像のうち直前に前記記憶領域に保存された撮影画像の前記特徴点と前記第１のステップにおいて前記記憶領域に保存された前記撮影画像の前記特徴点とをマッチングし、当該マッチングが失敗した場合には、さらに以前に前記記憶領域に保存された撮影画像の前記特徴点と前記第１のステップにおいて前記記憶領域に保存された前記撮影画像の前記特徴点とをマッチングする、請求項１に記載の顕微鏡画像情報処理方法。
前記撮影画像が前記記憶領域に保存されたことを検知することは、
ユーザから前記記憶領域に存在する監視対象フォルダの指定の入力を受け付け、
一定時間ごとに当該監視対象フォルダが更新されたか確認することで前記撮影画像が前記記憶領域に保存されたことを検知する、請求項１に記載の顕微鏡画像情報処理方法。
前記第４のステップにおいて生成される接合済み画像は、予め定められた比率で縮小された前記撮影画像を用いて生成される画像である、請求項１に記載の顕微鏡画像情報処理方法。
前記第３のステップにおいて、前記特定の条件下においては、前記記憶領域に保存されている全ての前記撮影画像のグローバル変換行列が算出され、前記特定の条件下ではない場合は、前記記憶領域に保存されている前記撮影画像のうちの特定の画像を基準とした、前記第１のステップにて保存された前記撮影画像のグローバル変換行列が算出される、請求項２に記載の顕微鏡画像情報処理方法。
前記特定の条件下において、前記グローバル変換行列はMax Spanning Treeによって算出される、請求項６に記載の顕微鏡画像情報処理方法。
前記特定の条件下において、前記第３のステップにおける前記マッチング処理が行われた全特徴点間の再投影誤差を最小二乗法を用いて最小化することによって前記グローバル変換行列を調整する、請求項７に記載の顕微鏡画像情報処理方法。
前記試料の撮影が終了するまで、前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像を含む領域について、当該領域と、全ての前記撮影画像と、をグローバル変換行列に従って配置するステップを、前記第５のステップの後にさらに含む、請求項１に記載の顕微鏡画像情報処理方法。
前記第６のステップにおいて、前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像に対してバンドル調整を実行してから、接合済み画像を複数の小領域に分割する、請求項１に記載の顕微鏡画像情報処理方法。
前記複数の小領域の各々について、前記第７のステップの後に、
前記第７のステップにおいて配置された前記処理対象の小領域と、前記特定された１または複数の前記撮影画像と、を接合する第８のステップと、
前記第８のステップにおいて接合された画像を一時ファイルとして記憶領域に保存する第９のステップと、が実行され、
前記複数の小領域の各々について前記第９のステップによって保存された複数の一時ファイルを記憶領域に順次書き込む第１０のステップと、
をさらに含む、請求項１に記載の顕微鏡画像情報処理方法。
前記バンドル調整は以下のレーベンバーグ・マーカート法で定義される式を評価することにより実行される、請求項１０に記載の顕微鏡画像情報処理方法。
ただし、
ｘ_ｉ：処理時点ｉにおける画像の変換行列の成分
Ｊ_ｉ：処理時点ｉにおけるＪａｃｏｂｉａｎ（右肩のＴは転置を示す）
ｒ_ｉ：処理時点ｉにおける特徴点間の誤差
λ：誤差の大きさに応じて調整されるゼロ以上の実数
Ｉ：単位行列
前記バンドル調整は以下のレーベンバーグ・マーカート法で定義される式を評価することにより実行される、請求項１０に記載の顕微鏡画像情報処理方法。
ただし、
ｘ_ｉ：処理時点ｉにおける画像の変換行列の成分および歪み補正パラメータ
Ｊ_ｉ：処理時点ｉにおけるＪａｃｏｂｉａｎ（右肩のＴは転置を示す）
ｒ_ｉ：処理時点ｉにおける特徴点間の誤差
λ：誤差の大きさに応じて調整されるゼロ以上の実数
Ｉ：単位行列
であり、前記歪み補正パラメータは、全処理時点における各画像に共通の歪みパラメータである。
前記試料は蛍光試料である、請求項１に記載の顕微鏡画像情報処理方法。
前記コンピュータシステムは、第１のコンピュータ装置と、前記第１のコンピュータ装置とは別体である第２のコンピュータ装置と、を含んで構成され、
前記第１のステップは前記第１のコンピュータ装置において実行され、
前記第２のステップ～前記第５のステップは前記第２のコンピュータ装置において実行される、請求項１に記載の顕微鏡画像情報処理方法。
顕微鏡を用いて観察される試料の一部分の撮影画像を取得して記憶領域に保存する第１のステップを実行し、
前記撮影画像が前記記憶領域に保存されたことを検知したことをトリガーとして、前記撮影画像の特徴点に関する情報である特徴点情報を算出する第２のステップを実行し、
前回までに前記記憶領域に保存された１または複数の前記撮影画像の前記特徴点情報と、前記第２のステップにおいて算出された前記撮影画像の前記特徴点情報と、を用いて、前回までに前記記憶領域に保存された１または複数の前記撮影画像の前記特徴点と前記第１のステップにおいて前記記憶領域に保存された前記撮影画像の前記特徴点とのマッチング処理を実行する第３のステップを実行し、
前記マッチング処理の結果と、前回までに前記記憶領域に保存された１または複数の前記撮影画像と、前記第１のステップにおいて前記記憶領域に保存された前記撮影画像と、を用いて接合処理を実行して接合済み画像を生成する第４のステップを実行し、
前記接合済み画像を表示出力する第５のステップを実行し、
前記試料を部分的に撮影するたびに前記第１のステップから前記第５のステップまでを繰り返し、
前記第１のステップから前記第５のステップが繰り返されることで生成された前記接合済み画像を含む領域を複数の小領域に分割する第６のステップと、
前記複数の小領域の各々について、
処理対象の小領域と交差する１または複数の前記撮影画像を特定し、
前記処理対象の小領域と、前記特定された１または複数の前記撮影画像と、をグローバル変換行列に従って配置する
ことを実行する第７のステップと、
をさらに実行する、顕微鏡画像情報処理システム。
請求項１～１５に記載の顕微鏡画像情報処理方法をコンピュータシステムに実行させる、コンピュータプログラム。