JP7099923B2 - 画像処理方法、コンピュータプログラムおよび記録媒体 - Google Patents

Description

この発明は、培養された細胞コロニーを撮像した画像から細胞コロニーを検出するための画像処理技術に関するものである。

例えば容器内で培養された細胞コロニーを含む試料に対し、コロニーの数や大きさなどの定量的な計測を行う目的で撮像が行われる場合がある。継続的な培養や観察等が必要とされる場合、非侵襲的な方法で撮像および計測を行うことが必要である。このような要求に応えるものとして、例えば特許文献１には、画像から心筋細胞のコロニーを抽出し、時間差をおいて撮像された複数の画像におけるコロニーの動きから拍動を検出する技術が記載されている。この技術においては、心筋細胞コロニーが有する形態的特徴量、例えばコロニーの輪郭や面積、形状、空間周波数、輝度等の情報を用いて、画像から細胞コロニーの領域が特定される。

特開２０１４－０７５９９９号公報

特に複数の細胞コロニーを含む画像においては、本来は独立して活動している複数の細胞コロニーが画像上において互いに重なりあるいは接触している場合があり得る。この場合において、上記した特徴量に基づく領域の抽出は必ずしも有効に機能しない場合がある。例えば、互いに接触した複数のコロニーからなるオブジェクトが１つのコロニーとみなされたり、複数のコロニーに囲まれた背景領域がコロニーの内部と誤認されたりするという誤判定が生じ得る。

このため、画像に複数の細胞コロニーを含み、しかもそれらが画像内で互いに接触しているような場合でも、それらを個別のコロニーとして確実に検出することのできるような画像処理技術の確立が求められる。

この発明は上記課題に鑑みなされたものであり、複数の画像コロニーを含んで撮像された画像から、各細胞コロニーを的確に分離して検出することのできる技術を提供することを目的とする。

この発明に係る画像処理方法の一の態様は、上記目的を達成するため、培養された細胞コロニーを複数含む試料を撮像した原画像を取得する工程と、前記原画像から、所定のエッジ強度を有する閉曲線を抽出する工程と、抽出された前記閉曲線で囲まれた領域を示す画像と前記原画像を２値化した画像との論理積に相当する合成画像を作成する工程と、前記合成画像に含まれるオブジェクトに対し分水嶺法によるオブジェクト分割を実行する工程とを備えている。

このように構成された発明では、原画像から所定のエッジ強度を有する閉曲線が抽出される。こうして抽出される閉曲線は試料中の細胞コロニーの輪郭を表している蓋然性が高いと言えるが、細胞コロニー同士が接触している場合、抽出される輪郭は複数の細胞コロニーを連結したものとなり本来の輪郭とは異なる。また、互いに接触する細胞コロニーで囲まれた背景領域も閉曲線の内部に取り込まれた形となり、この部分も細胞コロニーの内部と同じとみなされることとなってしまう。

一方、原画像を所定の濃度値（輝度値）を閾値として２値化した画像では、閾値を適宜に設定することにより、細胞コロニーの内部領域と外部つまり背景領域とを区分することが可能である。そのため、細胞コロニーで囲まれた背景領域と、その周囲の細胞コロニーの内部とを区別することができる。また、細胞コロニー内の内部テクスチャに関する情報をある程度維持することができる一方で、細胞コロニー内と背景との境界については閾値の設定によって大きく変動するため、輪郭の特定には向いていない。

そこで本発明では、原画像から抽出された閉曲線で囲まれる領域を表す画像と、原画像を２値化した画像との論理積に相当する合成画像を作成する。このようにすると、閉曲線抽出の結果および２値化の結果の両方において細胞コロニー内とみなされる領域のみが、細胞コロニーの内部領域として合成画像に現れることになる。これにより、連結した細胞コロニーで囲まれた背景領域が、背景として正しく取り扱われることとなる。

ただし、ここまでの処理では連結した細胞コロニーを個別のものとして互いに切り離すには至っていない。そこで、合成画像において細胞コロニーの内部とされるオブジェクトに対して分水嶺法によるオブジェクト分割がさらに実行される。これにより、画像内で連結する複数の細胞コロニーが分離される。

以上の処理により、画像内において複数の細胞コロニーが連結しているような場合であっても、それらは個別のものとして分離され、また連結した細胞コロニーで囲まれた背景領域についても正しく背景として細胞コロニーの内部とは分離することができる。これにより、本発明では、複数の細胞コロニーを含む原画像から、それらが接触している場合でも各細胞コロニーを分離して個別に検出することが可能である。

本発明に係る画像処理方法は、コンピュータをその実行主体とすることが可能なものである。この意味において、本発明は、コンピュータに上記処理を実行させるためのコンピュータプログラムとして実現することが可能である。また、当該コンピュータプログラムを記録した記録媒体として実現することも可能である。

上記のように、本発明の画像処理方法によれば、原画像から抽出された閉曲線に囲まれた領域を示す画像と原画像を２値化した画像との論理積を求め、さらに分水嶺法によるオブジェクト分割を実行することにより、画像内で複数の細胞コロニーが接触した状態であっても、それらを分離し個別のものとして検出することが可能である。

本実施形態の画像処理を示すフローチャートである。 １つのフレーム画像を模式的に例示する図である。 前処理工程を示すフローチャートである。 スフェア抽出工程を示すフローチャートである。 画像端補正工程を説明するための図である。 エッジ検出処理後の画像の例を示す図である。 輪郭領域画像の一例を示す図である。 ２値化画像の一例を示す図である。 論理積合成画像の一例を示す図である。 距離画像の一例を示す図である。 距離変換結果に基づく不要オブジェクト消去の概念を示す図である。 中心画像の一例を示す図である。 オブジェクト分割の概念を示す図である。 異なるフレーム間でのスフェアの対応関係を示す図である。 間引き処理の原理を示す図である。 動き検出工程を示すフローチャートである。 基準画像と比較画像との関係を示す図である。 基準画像と比較画像との差分を示す画像の例である。 探索範囲設定方法の第１の態様を示す図である。 探索範囲設定方法の第２の態様を示す図である。 収縮膨張画像と探索範囲画像との論理積合成画像を示す図である。 拍動判定工程を示すフローチャートである。 動き方向の判断例を示す図である。 この実施形態の画像処理を実行するコンピュータ装置の構成例である。

＜概要の説明＞
以下、本発明に係る画像処理方法の一実施形態について説明する。この実施形態における画像処理は、培養された心筋細胞のコロニーを含む試料を撮像した画像から細胞コロニーを検出し、さらには各細胞コロニーの拍動の有無を判定するための処理である。例えば適宜の培地中でｉＰＳ細胞（induced Pluripotent Stem Cell）やＥＳ細胞（Embryonic
Stem Cell）などの多能性幹細胞から分化培養された心筋細胞のコロニーを含む試料を、本実施形態の処理対象物とすることができる。

このような培養の初期段階では、幹細胞から分化した心筋細胞がいくつか集まって概ね球状のスフェアと呼ばれる小さなコロニーを形成する。試料には多数のスフェアが含まれ得るが、それらは互いに独立に拍動し、また中には拍動しないものもある。その他、未分化の細胞や細胞の代謝に伴い排出される老廃物等も試料には含まれ得る。以下に説明する本実施形態の画像処理は、このような画像から個々のスフェアを検出し、それらの拍動の有無を個別に判定することを目的とするものである。

この画像処理は、心筋細胞スフェアを含む試料が、所定の時間間隔で繰り返し撮像する動画撮像機能を備えた撮像装置によって撮像された画像に対してなされる。撮像装置が撮像後の画像データに対する後処理の１つとして本画像処理を実行してもよく、また撮像装置あるいは適宜の記憶手段から撮像画像に対応する画像データを受け取ったコンピュータ装置が、本画像処理を実行してもよい。

ここでは、汎用のコンピュータ装置が公知の撮像装置を制御して画像を取得し、該画像に対し本実施形態の画像処理を実行する態様を例示して説明する。撮像装置の構成については特に限定されず、培養された細胞等を含む試料を培地とともに動画撮像し、デジタル画像データとして出力する機能を有するものであればよい。例えば顕微鏡にデジタルビデオカメラを組み合わせた撮像装置を用いることができる。画像としては例えば明視野画像を好適に用いることができる。また、コンピュータ装置についても、例えばパーソナルコンピュータとして製品化されている一般的なハードウェア構成を有するものを利用可能である。

図１は本実施形態の画像処理を示すフローチャートである。まず、図１を参照してこの実施形態における画像処理の各工程の概略について説明し、その後で各工程における処理の具体的内容について詳しく説明する。なお、この画像処理は、撮像した画像からスフェアを検出するスフェア検出処理（ステップＳ１０１～Ｓ１０３）と、検出された各スフェアの拍動を検出する拍動検出処理（ステップＳ１０５～Ｓ１０７）との２つの処理に大別することができる。また、これらの間に、それぞれの処理に適したフレームレートを実現するための間引き処理が設けられる。

コンピュータ装置は、培地中で培養された心筋細胞のコロニーを含む試料を適宜のフレームレートで動画撮像する（ステップＳ１０１）。フレームレートとしては、心筋細胞の拍動周期として想定される周期よりは十分に短い撮像周期となることが好ましく、例えば一般的な動画撮像におけるフレームレートである３０ないし６０ＦＰＳ（Frame Per Second）程度とすることができる。また、撮像の継続時間としては、拍動する細胞コロニーの動きを１周期以上にわたり確実に捉えることのできる時間、例えば数秒とすることができる。

図２は１つのフレーム画像を模式的に例示する図である。フレーム画像Ｉｆでは、略一様かつ低濃度（高輝度）である背景ＢＧの中に、心筋細胞が略球状に集まって形成されたスフェアＳｐに対応するオブジェクトが背景より高い濃度を有して多数分布している。この他に、未分化細胞、遊離細胞、細胞から排出された老廃物の塊等、スフェア以外の物体Ｄｅに対応するオブジェクトも画像内には存在する。なお、背景とオブジェクトとの濃度の関係が上記と逆である場合には、以下の説明における濃淡についても逆に読み替えればよい。なお、図を見やすくするために、全てのオブジェクトに符号を付すことをせず、代表的な一部のオブジェクトにのみ符号を付している。以下の各図においても同様である。

図２に符号Ａ～Ｄを付して例示するように、複数のオブジェクトが画像内で互いに接し結合していることがある。このような結合が見られる理由としては、試料中で実際にオブジェクト同士が接触している場合と、試料内では離れているが画像内では接しているように見える場合とがある。

こうして撮像された画像については、フレームごとに、画像の微小な乱れを是正するための前処理工程が実行される（ステップＳ１０２）。そして、前処理後のフレーム画像に対してスフェア抽出工程が実行される（ステップＳ１０３）。この処理は、画像内に点在するオブジェクトの中からスフェアとみられるものを抽出し、その輪郭位置を特定するための処理である。これにより、スフェアが画像中でどの領域を占めているかが特定される。画像に複数のスフェアが含まれる場合、それらの輪郭が個別に特定される。

続いてフレーム間引きが行われる（ステップＳ１０４）。後述するように、本実施形態におけるスフェアの拍動検出は、フレーム画像間におけるスフェアの変位の大きさに基づいて行われる。ここで、試料の撮像はスフェアの拍動周期に対して十分に短い周期で実行されるが、撮像時刻の差が小さすぎる画像間ではスフェアの変位が小さくその検出が難しくなる。したがって全てのフレーム画像を用いてスフェアの動きを検出しようとすることは必ずしも現実的ではない。

一方で、撮像時刻の差が大きすぎると周期的な変位を見落とすおそれがある。拍動周期を考慮した適宜の間隔で画像を間引くことで、スフェアの変位をより確実に捕捉することが可能となる。また、画像処理のための演算量を大きく低減させることにも寄与する。この目的のために、撮像されたフレーム画像の中から所定のサンプリング間隔で一部を選択するフレーム間引きが行われる。

次に、画像間でのスフェアの動きが検出される（ステップＳ１０５）。上記したように、撮像時刻の異なる２つの画像間でのスフェアの変位から、当該スフェアの動きが検出される。そして、こうして検出された変位量から、各スフェアの拍動の有無が判定される（ステップＳ１０６）。画像内でのオブジェクトの変位は拍動によるものに限らない。検出された変位が拍動に起因するものであるか否かにより、拍動の有無が判定される。判定結果については適宜の形式で出力されユーザに提示される（ステップＳ１０７）。以上が本実施形態の画像処理の全体の流れである。以下、各処理の具体的内容につき順次説明する。

＜スフェア検出処理＞
図３は前処理工程を示すフローチャートである。撮像された各フレームの画像に対し、後述の画像処理の条件を一定とするためのコントラスト補正（ステップＳ２０１）およびノイズ除去（ステップＳ２０２）が行われる。これらの処理としては公知のものを適用可能であり、ここでは説明を省略する。例えばノイズ除去処理としてはバイラテラルフィルタ等の空間フィルタリング技術を用いることができる。なお、撮像条件によってはこれらの前処理を省くことができる場合がある。

次に、スフェア抽出工程（ステップＳ１０３）について説明する。スフェア抽出工程では、画像内に存在する各種オブジェクトの領域を特定し、それらの中からスフェアＳｐであるものを選択し、画像内で結合しているスフェアＳｐについては適切な境界で区分することで、個々のスフェアに対応する画像内のオブジェクトを特定する。この処理は撮像により得られたフレーム画像Ｉｆの全てについて実行される。

図４はスフェア抽出工程を示すフローチャートである。最初に１つのフレーム画像を選択し（ステップＳ３０１）、選択されたフレーム画像Ｉｆに対し画像端補正を行う（ステップＳ３０２）。図２に示すように、画像内のオブジェクトのうちフレーム画像Ｉｆの端部付近にあるものでは、その一部が画像端にかかり、全体が画像内に収まっていない場合がある。後述する処理においては輪郭とみなせる閉曲線で囲まれた領域を抽出する工程があるが、このように画像端で輪郭が途切れているオブジェクトの抽出漏れを回避するために、画像端補正が実行される。

図５は画像端補正工程を説明するための図である。図５（ａ）に示すように、フレーム画像Ｉｆの周縁部を取り囲むように、１画素分の幅を有する額縁状の領域Ｒｆを設ける。この領域の画素値（輝度値）は背景ＢＧと同程度、またはこれより低濃度（高輝度）となるように設定される。例えば濃度０（すなわち最大輝度の白色）とすることができる。この領域Ｒｆは、画像の端部に相当する画素の画素値を上記値に転換することによって設けられてもよく、また画像の端部の外側に上記画素値を有する仮想的な画素を付加することにより設けられてもよい。

処理後のフレーム画像Ｉｆでは、図５（ｂ）に示すように、画像端部の画素は全て背景濃度と同じまたはこれより低濃度となっている。このため、画像端にかかる（例えば図中に符号Ａ、Ｂで示す）オブジェクトも、画像内で閉じた輪郭を有することになる。これにより、輪郭の途切れに起因するオブジェクトの抽出漏れが回避される。なお図５（ｂ）では、画像端の様子を見やすくするためにフレーム画像Ｉｆの周縁部を点線で示している。

図４に戻って、こうして画像端が補正されたフレーム画像Ｉｆに対して、輝度値の変化が大きな箇所を検出するエッジ検出処理が実行される（ステップＳ３０３）。エッジ検出処理としては公知の方法、例えばSobelフィルタに代表される空間フィルタ処理を適用することができる。

図６はエッジ検出処理後の画像の例を示す図である。図２に示すフレーム画像Ｉｆに対応するエッジ検出処理後の画像（エッジ画像）Ｉｅでは、画像内のオブジェクトと背景との境界、つまりオブジェクトの外周部に相当する部分が強調されエッジとして検出される。図に符号Ａを付して例示するように、オブジェクトの内部テクスチャに起因するエッジも含まれ得る。また、図に符号Ｂを付して例示するように、予め画像端補正を行っているため、画像端にかかるオブジェクトについても閉じた輪郭が形成される。なお、検出された輪郭の微小な欠落を補償するために適宜の補正処理、例えば膨張収縮処理がさらに実行されてもよい。

こうして検出されたエッジが閉曲線をなすとき、該閉曲線は画像中のオブジェクトの輪郭である可能性が高い。そこで、エッジ画像Ｉｅにおいて閉曲線で囲まれた輪郭の内部領域を抽出しこれを塗りつぶすことで（ステップＳ３０４）、輪郭の内部と外部とを区分する。ここでは、塗りつぶされた輪郭内部の領域を「輪郭領域」と称し、塗りつぶし後の画像を「輪郭領域画像」と称することとする。

図７は輪郭領域画像の一例を示す図である。この輪郭領域画像Ｉｒは、図６のエッジ画像Ｉｅにおいて閉曲線で囲まれた領域を淡色（白）で、それ以外の領域を濃色（黒）で塗り分けた画像に相当する。なお、ここでは一般的な画像処理の例に倣いこのような塗り分けとしているが、濃淡が逆であってもよい。以下の処理においても同様である。

輪郭領域画像Ｉｒでは、閉曲線で囲まれたオブジェクトの内部が全て塗りつぶされるため、オブジェクトの内部テクスチャに対応する情報は失われる。また、図２に示す元のフレーム画像Ｉｆと対比すると、図７に符号Ａ、Ｂを付して例示するように、本来は背景部分であるが、周囲全てをオブジェクトによって取り囲まれているために塗りつぶされてしまう領域も生じる。すなわち、本来はオブジェクトの外部とされるべき背景領域がオブジェクトの内部として扱われることとなってしまう。

この問題を回避するための処理が次に実行される。すなわち、元のフレーム画像Ｉｆを適宜の閾値により２値化した画像が生成され（ステップＳ３０５）、さらにこの２値化画像と輪郭領域画像Ｉｒとを論理積合成した合成画像が生成される（ステップＳ３０６）。２値化されるフレーム画像Ｉｆとしては、前処理後のものが用いられることが好ましい。なお２値化のための閾値については固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としてもよい。

図８は２値化画像の一例を示す図である。より具体的には、図８は、図２のフレーム画像Ｉｆに前処理を施した画像を適宜の閾値で２値化し、閾値より高濃度の領域を淡色（白）で、より低濃度の領域を濃色（黒）で塗り分けることで得られる２値化画像Ｉｂを示している。２値化画像Ｉｂでは、オブジェクトの輪郭部分およびその内部の高濃度領域が白色で表される一方、元のフレーム画像Ｉｆにおいて低濃度であった、背景領域やオブジェクト内の一部の低濃度領域が黒色で表される。図８に符号Ａ、Ｂを付して例示するように、周囲をオブジェクトで囲まれた背景部分も黒色で表される。なお、画像Ｉｂの左上隅および右上隅の本来は背景である部分に現れている白色の領域は、背景濃度のムラが２値化の結果に影響を及ぼす可能性があることを象徴的に示したものである。

図９は論理積合成画像の一例を示す図である。図７に示す輪郭領域画像Ｉｒと図８に示す２値化画像Ｉｂとを論理積合成すると、図９に示す合成画像Ｉｓのように、元の両画像においてともに白色であった領域は白色で、それ以外の領域は黒色で表されることになる。これにより、図７の輪郭領域画像Ｉｒでは白色に塗りつぶされていた、オブジェクトで囲まれた背景領域が、他の背景部分と同じ黒色で表されることになる。一方で、２値化画像Ｉｂに含まれていたオブジェクトの内部テクスチャに起因する黒色領域は合成画像Ｉｓにおいても保存されている。

したがって、図９に示す合成画像Ｉｓは、オブジェクトに囲まれた背景を正しく背景として認識するという意味において、またオブジェクトの内部テクスチャに関する情報を保存するという意味において、図７の輪郭領域画像Ｉｒよりも有効なものであるといえる。また、合成画像Ｉｓは、明確なエッジとして検出された輪郭によって、背景の濃度ムラの影響を受けずオブジェクトと背景との境界が画定されるという意味において、図８の２値化画像Ｉｂよりも有効なものであるといえる。つまり、合成画像Ｉｓは、画像中のオブジェクトと背景とを区分するという目的において、輪郭抽出に基づく方法の持つ利点と、２値化に基づく方法の持つ利点とを併せ持ったものである。

ここまでの処理により、画像内でオブジェクトが占める領域と背景領域とが的確に区別される。しかしながら、前述したように、画像内のオブジェクトには、本実施形態の処理対象であるスフェア以外の不要なオブジェクトも含まれ、また複数のオブジェクトが画像内で接触、結合している場合がある。このため、不要オブジェクトの除去および結合したオブジェクトの分割のための処理が必要である。

この目的のために、画像内の各オブジェクトの構造解析を行う。具体的には、ステップＳ３０６で生成された合成画像Ｉｓに対し距離変換が実行される（ステップＳ３０７）。２次元２値画像における距離変換は、汎用の画像処理ソフトウェアでも実行可能な一般的な画像処理であり、画像内の各画素を当該画素から最も近い黒色画素までの距離によって表すというものである。

この実施形態において処理対象となる合成画像Ｉｓでは、背景に対応する画素において距離値０、背景画素に隣接しオブジェクトの輪郭に対応する画素において距離値１を取り、オブジェクト内部の画素ほど背景領域から遠ざかるため距離値は大きくなる。画像内の各画素の画素値を距離値に置き換え、例えば距離値が大きい画素ほど高い輝度値によって表すようにすれば、背景からの距離を濃淡で表した距離画像が得られる。

図１０は距離画像の一例を示す図であり、この距離画像Ｉｄは図９に示す合成画像Ｉｓを距離変換した結果を示すものである。距離変換により、合成画像Ｉｓは、背景からの距離が各画素の輝度値に反映された距離画像Ｉｄが得られる。距離画像Ｉｄは、背景部分が暗く、オブジェクトの周縁部から内部に向かうにつれて高輝度となるグレースケール画像である。したがって、面積が大きく円形に近いオブジェクトでは中央部の距離値が大きく、距離画像では輝度が高く（明るく）なる。一方、面積の小さいオブジェクトや細いオブジェクト、周囲の凹凸の大きいオブジェクト等では中央部の距離値は小さく、したがって輝度値が低く（暗く）なる。

処理対象物であるスフェアＳｐは、画像内で比較的大きく円形に近い形状を有する。したがって当該オブジェクトの内部には比較的大きな距離値を有する画素が存在する。一方、浮遊細胞等の小さなオブジェクトでは、オブジェクト内の各画素の距離値はより小さい。このことを利用して、画像内の各オブジェクトはスフェアＳｐとそれ以外の不要オブジェクトとに区別され、不要オブジェクトは距離画像Ｉｄおよび合成画像Ｉｓから消去される（ステップＳ３０８）。

図１１は距離変換結果に基づく不要オブジェクト消去の概念を示す図である。図１１（ａ）に示すように、距離画像Ｉｄの各画素の位置を座標（Ｘ，Ｙ）で表し、各画素の距離値を座標平面に対しプロットした場合を考える。なお、図ではＸ座標およびＹ座標を１つの座標軸で表しているが、より厳密にはＸＹ座標平面に対し距離値の軸が直交する形態となる。

このように表現したとき、距離値０の領域は背景に対応し、背景内で孤立した１つのオブジェクトは、０より大きい距離値を有する画素が連続する領域として表される。図１１（ａ）の例では、２つのオブジェクトＯｂｊ１、Ｏｂｊ２が互いに孤立して存在する。比較的サイズが大きく、また凹凸の少ないまとまった形状を有するオブジェクトＯｂｊ１では距離値のピークが大きく、より小さい細いオブジェクトＯｂｊ２ではピークが小さくなる。このことから、距離値に対し適宜の閾値Ｖｔｈを設定し、距離値のピーク値が閾値Ｖｔｈより大きいオブジェクトＯｂｊ１についてはスフェアとみなして画像に残す一方、距離値のピーク値が閾値Ｖｔｈより小さいオブジェクトＯｂｊ２については不要オブジェクトとみなして画像から消去する。

このようにすれば、比較的大きなオブジェクト、円形に近いオブジェクト等は画像に残る一方、小さなオブジェクトや細いオブジェクト、周囲の凹凸の激しいオブジェクト等については画像から消去される。図１１（ｂ）は不要オブジェクト消去後の距離画像Ｉｄａの一例を示している。符号Ａ～Ｄを付した位置では、図１０に現れていた小さなオブジェクトや細長いオブジェクトが現れていない。すなわち、このようなオブジェクトは消去されている。画像の状態に応じた的確な不要オブジェクト除去を可能とするために、上記閾値Ｖｔｈについては固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としておくことが好ましい。

このようにして、不要オブジェクトを除外した距離画像とこれに対応する中心画像Ｉｃａとを作成することができる（ステップＳ３０９）。ここでいう中心画像Ｉｃａは、距離画像Ｉｄａにおいて、距離値が極大となり、かつその極大値が所定値より大きい画素の位置を表す画像である。

図１２は中心画像の一例を示す図である。より具体的には、図１２（ａ）は図１１（ｂ）に示す距離画像Ｉｄａに対応する中心画像Ｉｃａである。また、図１２（ｂ）は、図１２（ａ）の中心画像Ｉｃａを、図９に示す合成画像Ｉｓに重畳した参考図である。なお、図１２（ｂ）においては、中心画像Ｉｃａにより示される各点の視認性を高めるため、その位置を黒丸印で示す。

これらの図からわかるように、また図１２（ｂ）に符号Ａ，Ｂを付して例示したように、円形に近いオブジェクトにあっては、このような画素は当該オブジェクトの中心に近い位置にある。この意味において、中心画像は各オブジェクトの中心位置を表しているといえる。そこで、中心画像Ｉｃａに現れる点を「中心点」と称することとする。図１２（ｂ）の画像には中心点を含まないオブジェクトも含まれているが、これらは不要オブジェクトとみなされるものであり、中心点の導出対象から除外されている。

図１２（ｂ）に符号Ｃ，Ｄを付して例示するように、複数のオブジェクトが結合しているとみられる複雑な形状のオブジェクトや、符号Ｅを付した顕著なくびれを有するオブジェクトでは、１つのオブジェクトに複数の中心点が存在する。言い換えれば、このように複数の中心点が存在するオブジェクトは、中心点の数に対応する複数のオブジェクトが結合したものと考えられる。

そこで、複数の中心点が存在するオブジェクトについては、それらの中心点の間に適宜の境界線を設定することで、それぞれ中心点を１つずつ有するオブジェクトに分割する。こうすることで、結合したオブジェクトを分割することができる。このような考え方は、「分水嶺（Watershed）法」と呼ばれるオブジェクト分割のための画像処理のものに他ならない。すなわち、距離変換と不要オブジェクト消去とにより得られた距離画像Ｉｄａおよび中心画像Ｉｃａの情報に基づき、合成画像Ｉｓに対し分水嶺法によるオブジェクト分割を実行することで（ステップＳ３１０）、画像内で結合するオブジェクトが分割される。

図１３はオブジェクト分割の概念を示す図である。図１３（ａ）に模式的に示すように、２つの中心点Ｃ１、Ｃ２を有するオブジェクトＯｂｊ３は、それらの中心点Ｃ１、Ｃ２の間に、破線で例示する境界線を設定することで２つに分割される。この境界線を適切に設定するための画像処理アルゴリズムの１つが分水嶺法である。３以上の中心点を含むオブジェクトについても同様にすることができる。

図１３（ｂ）はオブジェクト分割後の画像の例であり、より具体的には、図１１（ｂ）の距離画像Ｉｄａおよび図１２（ａ）の中心画像Ｉｃａの情報に基づく分水嶺法により、図９に示す合成画像Ｉｓをオブジェクト分割して得られたオブジェクトの輪郭を示す画像Ｉｏである。ただし、予め合成画像Ｉｓから不要オブジェクトは消去されているものとする。符号Ａ～Ｃを付して例示するように、合成画像Ｉｓでは結合していた複数のオブジェクトが分離され、また符号Ｄを付して例示するように、顕著なくびれを有するオブジェクトもそのくびれ部分から２つに分割される。こうして得られた１つ１つのオブジェクトは、それぞれ１つのスフェアとみなすことができる。

このようにして１つのフレーム画像Ｉｆから抽出されたスフェアＳｐについてはインデックスが作成される（ステップＳ３１１）。例えば各スフェアＳｐに付した固有の識別記号と、その代表点（例えば中心）の位置を表す情報とを関連付けたインデックスを作成することができる。その他にも、例えば各スフェアＳｐの面積や輪郭の形状、周囲長などを表す情報をインデックスに含めてもよい。

こうして１つのフレーム画像Ｉｆについてスフェア抽出が終了する。次に処理すべきフレーム画像がある場合には（ステップＳ３１２においてＹＥＳ）、ステップＳ３０１からの処理を繰り返すことで、各フレーム画像が順次処理される。最後のフレーム画像まで処理が終了すれば（ステップＳ３１２においてＮＯ）、フレーム抽出処理は完了する。

上記したステップＳ３０１～Ｓ３１１の処理を各フレーム画像Ｉｆについて一律に実行することで、当該フレーム画像Ｉｆ中のスフェアＳｐを抽出することが可能である。この場合、複数のフレーム画像からそれぞれ抽出されるスフェアの間に何ら関連付けはなされない。しかしながら、実際には、小さな時間差で撮像された２つの画像でそれぞれ抽出されるスフェアは、ほぼ１対１の関係で関連付けることができる。言い換えれば、１つの画像に含まれるスフェアの各々は、多少の位置や形状の変化はあるとしても、他の１つの画像にもほぼ同様の位置に同じような形状で存在しているはずである。

スフェアの拍動を検出するという目的に鑑みれば、１つのスフェアの位置が経時的にどう変化するかを追跡する必要がある。したがって、１つの画像に含まれるスフェアが他の画像中のどのスフェアに対応するかを把握しなければならない。このような観点から、第２フレーム以降の画像に対する処理は以下のように改変することができる。

すなわち、改変されたスフェア抽出処理において、第２フレームの画像に対する処理では、ステップＳ３０９として記載された距離画像および中心画像の生成が省略される。そして、ステップＳ３１０の分水嶺法によるオブジェクト分割においては、第１フレームの画像に対する処理で生成された距離画像および中心画像の情報が使用される。つまり、第２フレームの画像中の各オブジェクトについては、その中心点が第１フレーム画像中のオブジェクトと同じ位置にあるものとみなされる。このようにする理由につき以下に説明する。

図１４は異なるフレーム間でのスフェアの対応関係を示す図である。図１４（ａ）は画像中で孤立したスフェアの例を示している。なお、図１４（ａ）はフレーム画像の一部のみを示しており、また点線による格子は画像内でのスフェアの位置を把握しやすくするために加えた仮想的なものである。図左に示す第Ｎフレームの画像に含まれるスフェアＳｐは、図右に示す第（Ｎ＋１）フレームの画像では少し変位した位置に現れる。このため、第（Ｎ＋１）フレーム画像におけるスフェアＳｐの中心点Ｃ（Ｎ＋１）は、第Ｎフレーム画像におけるスフェアＳｐの中心点Ｃ（Ｎ）の位置から少し移動していることがある。

しかしながら、短い撮像間隔の間での変位量は限定的であり、第Ｎフレーム画像における中心点Ｃ（Ｎ）は、第（Ｎ＋１）フレーム画像におけるスフェアＳｐの内部に含まれている、またはそれに極めて近い位置にあると予想される。したがって、第（Ｎ＋１）フレーム画像において中心点Ｃ（Ｎ）を含むスフェアが、または第（Ｎ＋１）フレーム画像において中心点Ｃ（Ｎ）に最も近い位置にあるスフェアが、第Ｎフレーム画像におけるスフェアと同一であるといえる。このことから、異なるフレーム間で対応するスフェアを相互に関連付けることができる。

図１４（ｂ）は画像中で結合したスフェアの例を示している。第Ｎフレームの画像に２つのスフェアＳｐ１，Ｓｐ２が含まれ、これらが画像内で結合しているとする。スフェアＳｐ１，Ｓｐ２それぞれの中心点をＣ１（Ｎ），Ｃ２（Ｎ）とする。また、このときにオブジェクト分割により画定されるスフェア境界線を破線Ａによって表す。

第（Ｎ＋１）フレームの画像において、図に点線で示すように、スフェアＳｐ２が右方向に変位したとする。このときのスフェアＳｐ２の中心点をＣ２（Ｎ＋１）とすると、スフェアＳｐ２の動きに対応して中心点Ｃ２（Ｎ＋１）も右方向に変位していると考えられる。一方、スフェアＳｐ１はほとんど動かず、第（Ｎ＋１）フレームにおける中心点Ｃ１（Ｎ＋１）は、第Ｎフレームにおける中心点Ｃ１（Ｎ）とほぼ同じである。

第（Ｎ＋１）フレームの画像において、検出される中心点Ｃ１（Ｎ＋１），Ｃ２（Ｎ＋１）の位置に基づきオブジェクト分割を行ったとき、破線Ｂで示すように、中心点間の距離の増大に起因して境界線も右方向にシフトすることになる。これにより、スフェアＳｐ１が右方向に広がったように見える。しかしながら、実際にはスフェアＳｐ２が動いたのであって、スフェアＳｐ１が占める領域はほぼ変わらない。このようなオブジェクト分割の結果と実態との乖離を抑えるには、第（Ｎ＋１）フレームにおいても、第Ｎフレームで特定された中心点Ｃ１（Ｎ），Ｃ２（Ｎ）の情報を用いるようにすればよい。

第（Ｎ＋２）フレーム以降の画像についてオブジェクト分割を行う際に、どのフレームで特定された中心画像の情報を用いるかについてはいくつかの考え方があり得る。例えば、全てのフレームの処理に対して単一のフレーム、例えば第１フレームにおける中心画像の情報を適用することが考えられる。また、第（Ｎ＋１）フレームの処理に対してその直前に撮像された第Ｎフレームにおける中心画像の情報を適用することが考えられる。さらに、一定のフレーム数ごとに中心画像を更新し、次の更新までは同じ情報を用いて各フレームの処理を行うようにしてもよい。

このように、一のフレーム画像に対するオブジェクト分割に際し、異なる時刻に撮像されたフレーム画像の距離情報および中心情報を用いて分水嶺法を実行することにより、それらのフレーム画像間において対応するスフェアの関連付けを行うことができる。例えば各フレーム画像で検出されたスフェアにつきインデックスを作成する際に、対応するスフェアに共通の識別記号を付すようにすれば、画像間でのスフェアの対応関係を明確に記録することができる。

以上のように、本実施形態のスフェア検出処理では、連続撮像された各フレーム画像からスフェアを検出するのに際して、次のような要素技術を含んでいる。まず、画像端にかかるオブジェクトが除外されるのを防止するために、画像の周縁部を背景とみなせる濃度に設定し、画像端にかかるオブジェクトが閉じた輪郭で囲まれるようにする。

また、エッジ検出により元のフレーム画像から抽出されたオブジェクトの輪郭に基づく輪郭領域画像と、フレーム画像を濃度（輝度）に基づいて２値化した２値化画像との論理積合成画像が、画像においてスフェア等のオブジェクトと背景との境界を示すものとされる。これにより、オブジェクトに取り囲まれた背景領域が正しく背景として認識されるようになる。

そして、距離変換の結果に基づきスフェアでない可能性の高いオブジェクトを消去し、残るオブジェクトに対し分水嶺法によるオブジェクト分割を実行することで、画像内で結合したスフェアをそれぞれ個別のオブジェクトとして特定することができる。このとき、撮像時刻の異なる画像における距離変換の結果を用いてオブジェクト分割を行うことで、経時的に動きのあるスフェアを画像間で関連付けることができ、また動きのないスフェアと動きのあるスフェアとが結合している場合でも、それらを的確にかつ安定的に分割することができる。

＜間引き処理＞
図１５は間引き処理の原理を示す図である。図１５（ａ）に示すように、試料を所定のフレームレートで連続撮像することにより、一定の時間間隔で撮像された複数の原画像Ｉｇ１，Ｉｇ２，…が取得される。例えばフレームレートが６０ＦＰＳであるとき、時系列において隣り合う、つまり撮像時間間隔が最も小さい２つの原画像間の撮像時間間隔は（１／６０）秒である。一方で、心筋細胞スフェアの拍動周期は概ね１秒程度であるため、隣り合う２つの原画像の間では拍動に起因するスフェアの変位は極めて小さい。拍動に起因するスフェアの変位を効率よくかつ確実に検出するためには、拍動周期を考慮して画像の間引きを行うことが有効である。

例えば１秒の周期で振動するスフェアの変位を検出するためには、フレーム間の撮像時間間隔は１／６秒ないし１／３秒程度が好ましい。この場合、原画像を１／１０ないし１／２０程度に間引くことになる。この間引きの度合いを以下では「間引き率」と称する。種々の試料に柔軟に対応するために、間引き率については固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としておくことが好ましい。なお、この場合の原画像Ｉｇ１，Ｉｇ２，…としては、例えばスフェア検出工程における前処理後のフレーム画像を用いることができる。

図１５（ｂ）に示すように、原画像Ｉｇ１，Ｉｇ２，…から、間引き率に応じた一定のサンプリング間隔でサンプリング抽出した画像を時系列に沿って再配置することで、間引きが実現される。間引き後の画像群Ｉｆ１，Ｉｆ２，Ｉｆ３，…は、元のフレームレートに間引き率の逆数を乗じた値を実効的なフレームレートとする連続画像とみなすことができる。以下の処理は、これらの各画像を実質的なフレーム画像として実行される。そこで、以下の処理の説明において単に「フレームレート」および「フレーム画像」というとき、これらは間引き後のものを指すこととする。

＜拍動検出処理＞
次に、上記のようにして検出された各スフェアについての拍動検出処理について説明する。本実施形態の拍動検出処理は、図１に示すように、動き検出工程（ステップＳ１０５）および拍動判定工程（ステップＳ１０６）を備えている。これらの工程は、画像中のスフェアの各々について実行される。

図１６は動き検出工程を示すフローチャートである。また、図１７は基準画像と比較画像との関係を示す図である。この動き検出工程においては、間引き後のフレーム画像Ｉｆ１，Ｉｆ２，…から選択した１つの画像を基準画像とし（ステップＳ４０１）、これと撮像時刻が近いが異なる他の１つの画像を比較画像として（ステップＳ４０２）、両画像間におけるスフェアの変位を求め、その結果に基づき両画像の撮像時刻の間でのスフェアの動きを検出する。まず、基準画像と比較画像との関係につき、図１７を参照して説明する。

スフェアの拍動周期に応じて適宜に間引かれたフレーム画像Ｉｆ１，Ｉｆ２，…における動き検出では、まず最初の画像Ｉｆ１が基準画像とされる。そして、これより後に撮像された画像Ｉｆ２，Ｉｆ３，Ｉｆ４を順次比較画像に設定し、それぞれ基準画像Ｉｆ１との間でスフェアの変位が求められる。

時系列において隣り合う画像間での比較により、フレーム画像間の撮像時間間隔に等しい時間におけるスフェアの動きがわかる。ただし、スフェアの拍動周期は必ずしも一定ではなく、またスフェアごとに独立である。このため、基準画像と比較画像との撮像時間の差が、拍動に起因するスフェアの動きを捕捉するのに最適なものでない場合があり得る。例えば、拍動するスフェアが２つの画像の間で偶然に同じ位置にあるということが起こり得る。このことに起因する検出漏れを回避するために、比較画像を順次入れ替え、数種類の時間間隔でスフェアの変位を検出する。そして、基準画像と比較画像とのいずれかの組み合わせにおいてスフェアの変位が認められると動きありと判断するようにすれば、より確実な動き検出が可能となる。

ここでは１つの基準画像に対し３つの比較画像を設定することとしているが、比較画像として用いる画像の数はこれに限定されず任意である。ただし、比較画像の数を増やすと演算処理量も増大し実行装置の負荷が大きくなる、そのため、試料の実状に応じて比較画像の数と間引き率とを適切に設定することで、上記のような検出漏れの防止を図ることが望まれる。

基準画像Ｉｆ１に対し３つの比較画像を適用した後、次に基準画像を画像Ｉｆ２に切り替え、これに対しても同様に３つの画像Ｉｆ３，Ｉｆ４，Ｉｆ５を順次比較画像に設定してスフェアの変位が求められる。その後も同様に、基準画像を画像Ｉｆ３，Ｉｆ４，…と順番に切り替えながら、その都度３つの比較画像を順次設定して変位を求める。図１６のフローチャートでは、ループ処理によってこれが実現されている。このようにすることで、種々の周期で拍動するスフェアの動きを漏れなく検出することが可能となる。

図１６に戻り動き検出工程の処理内容についての説明を続ける。一の基準画像および一の比較画像が選出されると（ステップＳ４０１，Ｓ４０２）、両画像間のコントラストの違いに起因する検出誤差を低くするために、コントラスト補正が実行される（ステップＳ４０３）。補正処理としては公知のものを適用可能である。また、補正が必要なければ省くことも可能である。

続いて、基準画像と比較画像との差分画像が生成される（ステップＳ４０４）。差分画像はさらに２値化される（ステップＳ４０５）。両画像間でスフェアに動きがなければ画像内容に変化はなく、差分はゼロとなるはずである。一方、変位が大きければ差分も大きくなる。したがって、差分画像を２値化した画像は変位の大きさを定量的に表すための情報を有している。

図１８は基準画像と比較画像との差分を示す画像の例である。互いに対応するスフェアＳｐを含む基準画像Ｉ１と比較画像Ｉ２とを考える。図において点線による格子は画像内でのスフェアの位置を把握しやすくするために加えた仮想的なものである。基準画像Ｉ１と比べると、比較画像Ｉ２におけるスフェアＳｐは、その外形形状を変化させながらわずかに右方向に変位している。これらの差分を求めると画像Ｉ３が得られる。差分画像Ｉ３では、スフェアＳｐの外周部の変位に起因する概ね環状の像と、その内部に含まれるスフェア内部のテクスチャの変化に起因する像とが現れている。実際の差分画像は多階調画像であり、定量的な取り扱いを容易にするために２値化することで、２値化画像Ｉ４が得られる。

２値化画像Ｉ４に対し、公知の収縮膨張処理が実行される（ステップＳ４０６）。これにより、２値化画像Ｉ４から微細な斑点状のノイズ成分が除去される。なお、以下ではノイズ除去後の２値化画像を「収縮膨張画像」と称し符号Ｉ４ａにより表すこととする。この２値化画像Ｉ４ａは、スフェアの外周部の変位に起因する差分像とスフェア内部のテクスチャ変化に起因する差分像とを含む。このうちテクスチャ変化に起因する像は、スフェアの変位量を定量的に求めるに当たっては誤差要因となるものである。その理由は以下の通りである。

内部テクスチャの変化は、変位に直接関係のない細胞の活動によっても生じる。このため、スフェアが変位せず内部テクスチャのみが変化する場合がある。このような変化が変位と誤認されることがないようにする必要がある。また、仮に内部テクスチャに変化がないとしても、スフェア全体が変位することによって、同一位置の画素同士の差分を求めれば有意な差が現れることとなる。この差の大きさはスフェアの内部構造に依存するから、差分の値は変位量のみならずスフェア内部構造にも影響を受けていることになる。より具体的には、同量の変位であってもスフェアの内部構造によって差分の大きさが変わる。

したがって、スフェアの変位のみを定量的に求めるためには、スフェアの内部構造の影響を排除し、スフェア外周部の変位に起因する差分のみを抽出する必要がある。このために、本実施形態では、スフェア外周部の近傍に探索範囲を設定し（ステップＳ４０７）、この探索範囲内での画像間の差分を有意な値として求めるという手法が採用されている。以下、探索範囲設定方法の２つの態様につき説明するが、いずれの方法が適用されてもよい。

図１９は探索範囲設定方法の第１の態様を示す図である。より具体的には、図１９（ａ）は探索範囲設定方法の第１の態様を示すフローチャートであり、図１９（ｂ）は処理過程における画像の例を模式的に示す図である。図１９（ａ）に示すように、基準画像中のスフェアおよび比較画像中のスフェアのそれぞれにつき輪郭描画が行われる（ステップＳ５０１、Ｓ５０２）。各フレーム画像におけるスフェアの輪郭は、先のスフェア抽出処理で既に求められている。こうして求められた輪郭情報に基づき輪郭描画が行われるが、その線幅は１画素より大きな値に設定される。つまり、輪郭を膨張させた状態で描画する。

輪郭描画の結果、図１９（ｂ）に示すように、基準画像および比較画像それぞれにおいて特定されたスフェアの輪郭（点線）を所定の線幅に膨張させて描画した輪郭画像Ｉ５、Ｉ６が得られる。これらを論理和合成すると（ステップＳ５０３）、図１９（ｂ）下部に示すように、合成された画像Ｉ７には、基準画像および比較画像それぞれにおけるスフェアの輪郭を共に含み、かつそれらの輪郭の周囲まで広がる環状の帯状領域が現れる。この帯状領域が探索範囲Ｒｓとされる。ここでは画像Ｉ７を第１の態様における「探索範囲画像」と称することとする。

外周部の位置の変化を確実に差分として捕捉するためには、膨張させる輪郭の線幅は太い方が望ましい。その一方で、そうして探索範囲を広くすると、スフェアの外周部の変化以外の変化、例えば外周部に近い領域の内部テクスチャの変化も検出されてしまうことになる。そこで、上記線幅については固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としておくことが好ましい。

図２０は探索範囲設定方法の第２の態様を示す図である。より具体的には、図２０（ａ）は探索範囲設定方法の第２の態様を示すフローチャートであり、図２０（ｂ）は処理過程における画像の例を模式的に示す図である。なお、概念の理解を容易にするために、図２０（ｂ）ではスフェアの外形として単純な楕円形が用いられている。この方法では、先に求められている輪郭情報に基づき、基準画像、比較画像のそれぞれにつき、スフェアの輪郭内部を表す画像Ｉ８，Ｉ９を描画する（ステップＳ６０１、Ｓ６０２）。そして、それらの画像Ｉ８，Ｉ９は排他的論理和合成される（ステップＳ６０３）。

これにより得られる画像Ｉ１０では、２つの画像の一方でのみスフェアの輪郭内とされる領域のみが淡色（白）で表され、それ以外の領域は濃色（黒）で表される。この淡色領域を所定幅だけ膨張させる膨張処理を実行すると（ステップＳ６０４）、処理後の画像Ｉ１１には、基準画像および比較画像それぞれにおけるスフェアの輪郭（点線）を共に含み、かつそれらの輪郭の周囲まで広がる環状の帯状領域が現れる。この帯状領域が探索範囲Ｒｓとされる。ここでは画像Ｉ１１を第２の態様における「探索範囲画像」と称することとする。

こうしていずれかの態様の探索範囲設定方法により生成される、探索範囲Ｒｓを表す探索範囲画像をマスクとして、収縮膨張画像Ｉ４ａに作用させる。このことは、収縮膨張画像Ｉ４ａと探索範囲画像Ｉ７または探索範囲画像Ｉ１１とを論理積合成することと等価である（ステップＳ４０８）。

図２１は収縮膨張画像と探索範囲画像との論理積合成画像を示す図である。ここでは探索範囲画像として第１の態様の画像Ｉ７を用いて説明する。基準画像と比較画像との間におけるスフェアの外周部および内部の差分を表す収縮膨張画像Ｉ４ａと、探索範囲画像Ｉとを論理積合成することで、合成画像Ｉ１２では、スフェアの外周部の近傍における変化に対応する像のみが示され、収縮膨張画像Ｉ４ａに存在していたスフェア内部の像はマスクされて合成画像には現れない。これにより、外周部の位置変化、つまりスフェアの変位のみを選択的に抽出することが可能になる。この合成画像Ｉ１２を以下では「動き差分画像」と称する。

変位の定量化は以下のようにして行われる。すなわち、動き差分画像Ｉ１２から、動き差分面積が求められる（ステップＳ４０９）。動き差分面積は、動き差分画像Ｉ１２における白色部分の面積として求められ、簡易的には該画像における白色画素の数により表すことができる。基準画像と比較画像との間におけるスフェアの変位が大きいほど動き差分面積は大きくなり、したがって動き差分面積は、変位の大きさを定量的に表す数値として利用可能である。

ただし、動き差分面積はスフェアの大きさの影響を受ける。すなわち、面積が大きく外周部の周囲長の長いスフェアとより小さいスフェアとでは、たとえ変位量が同じであったとしても、変位の結果として現れる動き差分面積は、大きいスフェアの方が当然に大きくなる。これを是正し、スフェアの大きさに関わらず変位の大きさを表すことのできる指標として、以下の式：
Ｓ＝（動き差分面積）／（スフェアの平均輪郭長さ）
により、動きスコアＳを定義する。ここに、「スフェアの平均輪郭長さ」は、基準画像におけるスフェアの輪郭長さと比較画像におけるスフェアの輪郭長さとの平均値である。スフェアの輪郭長さについては、スフェア抽出処理の結果から導出することができ、例えば輪郭に相当する画素の数により表すことができる。

このように定義される動きスコアＳは、スフェアの変位により画像間に生じる差分がスフェアのサイズで正規化されたものであり、スフェアの大小に関わりなくその変位の大きさを指標するものとなる。

基準画像と比較画像との一の組み合わせごとに、画像内の全てのスフェアにつき動きスコアＳが算出される（ステップＳ４１０）。そして、上記したように、基準画像と比較画像との組み合わせを順次変化させながら（ステップＳ４１１、Ｓ４１２）、それぞれの組み合わせでスフェアの動きスコアＳが算出される。このうち動きスコアＳが予め定められた閾値より大きいものについては、当該スフェアが有意な動きを示したものとみなすことができる。動きの検出されたスフェアについては、その動きの方向とともに記録されることが望ましい。

図１に戻って、次に、動きスコアＳが算出された各スフェアにつき、その算出結果に基づき拍動判定が行われる（ステップＳ１０６）。スフェアの拍動の有無については、例えば以下のようにして判定することができる。

図２２は拍動判定工程を示すフローチャートである。より具体的には、図２２は１つのスフェアについて拍動の有無を判定するための処理を示している。各スフェアには固有の識別番号を含むインデックスが付されているので、各フレームでの検出結果から同一識別番号のものを抽出することにより、１つのスフェアの動きを把握することができる。動き検出工程で求められた１つのスフェアについての動き検出結果を時系列順に配列すると（ステップＳ７０１）、当該スフェアが経時的にどのような動きを示したかが明らかになる。

拍動はスフェアの継続的かつ周期的な運動である。したがって、スフェアが拍動していれば、当該スフェアについては複数回の動きが検出され、かつそれらの中には変位方向が反対のものが含まれているはずである。そこで、判定対象のスフェアがこのような条件を満たしていれば（ステップＳ７０２およびステップＳ７０３においてともにＹＥＳ）、当該スフェアは拍動していると判定される（ステップＳ７０４）。一方、動きの検出が１回のみである（ステップＳ７０２においてＮＯ）、あるいは複数回の動きの方向が同じである（ステップＳ７０３においてＮＯ）場合には、当該スフェアは拍動していないと判定される（ステップＳ７０５）。このようにして、各スフェアの拍動の有無を判定することができる。動き方向が反対であるか否かについては、例えば以下のようにして判断することができる。

図２３は動き方向の判断例を示す図である。検出された動きの方向をベクトル表示したとき、図２３（ａ）に示すように、１つのベクトルＶ１の方向と、他の１つのベクトルＶ２の方向とのなす角θが９０度より小さければ、これらの動きは同方向のものと判断することができる。このような動きは拍動によるものではなく、単にスフェアが一方向へ移動しているものと考えられる。一方、図２３（ｂ）に示すように、２つのベクトルＶ１，Ｖ２の角θが９０度より大きいとき、これらの動きは反対方向のものと判断することができる。なおスフェアの動き方向については、上記した動き検出結果から求めてもよく、また例えばスフェアの重心位置の経時的な変化から求めてもよい。

各スフェアの拍動に関する判定結果は、任意の態様で出力される（ステップＳ１０７）。例えば図１３（ｂ）に示すスフェア輪郭の画像Ｉｏにおいて、各スフェアの拍動の有無を識別できるような加工、例えば色分けを施して表示出力することが考えられる。また、図２に示すフレーム画像Ｉｆに対して、拍動するスフェア、拍動しないスフェアおよびスフェア以外のオブジェクトを互いに区別する加工を施してもよい。また、各スフェアの拍動の有無やその振幅、周期などの詳細な情報をインデックス情報とともにリストとして出力する態様であってもよい。その他、用途に応じて種々の態様で判定結果を出力することが可能である。

以上のように、この実施形態の拍動検出処理では、撮像時刻が異なる基準画像と比較画像との画像内容の差分に基づいてスフェアの変位が求められ、その結果からスフェアの動きが検出される。拍動に起因する変位はスフェアの外周部の近傍に現れる一方、スフェア内部のテクスチャ変化は必ずしも拍動の有無を反映しない。このため、本実施形態では、抽出されたスフェアの輪郭の近傍に探索範囲を設定し、該範囲内での画像の差分を変位の大きさの指標としている。このため、内部テクスチャの変化に影響されることなく、拍動に起因するスフェアの変位を効率よく求めることが可能である。

また、スフェアの拍動周期が必ずしも一定でなく、個体ばらつきもあることから、基準画像に組み合わせる比較画像については撮像時刻の異なる複数の画像を順次適用している。このため、種々の周期で生じ得るスフェアの拍動を確実に検出することが可能である。

またスフェアの変位の大きさについては、変位によって差分画像に現れる領域の面積をスフェアの輪郭長さで正規化した動きスコアＳを用いて評価する。このため、スフェアの大きさの違いに対してロバストな検出および判定が可能となっている。

＜装置構成例＞
図２４はこの実施形態の画像処理を実行可能なコンピュータ装置の構成例である。コンピュータ装置１は、例えばパーソナルコンピュータとして一般的な構成を有するものであり、ＣＰＵ（Central Processing Unit）１０、メモリ１４、ストレージ１５、入力デバイス１６、表示部１７、インターフェース１８およびディスクドライブ１９などを備えている。

ＣＰＵ１０は、予め用意された制御プログラムを実行することで、上記した画像処理を実行するための機能ブロックとしての画像処理部１１をソフトウェア的に実現する。なお、画像処理部１１を実現するための専用ハードウェアが設けられてもよい。メモリ１４はＣＰＵ１０の演算過程で生成される各種データを一時的に記憶する。ストレージ１５は、ＣＰＵ１０が実行すべき制御プログラムのほか、原画像の画像データや処理後の画像データ等を長期的に記憶する。

入力デバイス１６は、オペレータからの指示入力を受け付けるためのものであり、例えばマウス、キーボードなどを含む。また、表示部１７は画像を表示する機能を有する例えば液晶ディスプレイであり、原画像や処理後の画像、オペレータへのメッセージ等種々の情報を表示する。なお、入力デバイスと表示部とが一体化されたタッチパネルが設けられてもよい。

インターフェース１８は、電気通信回線を介して外部装置との間で各種データ交換を行う。ディスクドライブ１９は、画像データや制御プログラム等各種のデータを記録した外部の記録ディスク２を受け入れる。記録ディスク２に記憶された画像データや制御プログラム等は、ディスクドライブ１９により読み出され、ストレージ１６に記憶される。ディスクドライブ１９はコンピュータ装置１内で生成されたデータを記録ディスク２に書き込む機能を備えていてもよい。

本実施形態の画像処理をコンピュータ装置１に実行させるための制御プログラムについては、これを記録した記録ディスク２にディスクドライブ１９がアクセスして読み出される態様であってもよく、インターフェース１８を介して外部装置から与えられる態様であってもよい。原画像データについても同様である。

＜その他＞
なお、本発明は上記した実施形態に限定されるものではなく、その趣旨を逸脱しない限りにおいて上述したもの以外に種々の変更を行うことが可能である。例えば、上記実施形態では、試料の撮像、スフェアの抽出および拍動の判定までが一連の処理として記載されているが、常にこのような処理を要するものではない。例えば、過去に撮像された試料の画像を用いてスフェア抽出および拍動判定を行うようにしてもよい。また、スフェア抽出と拍動判定とについてもそれぞれ個別に実行されてもよい。

また、上記実施形態では、画像内におけるスフェアの変位を求めることにより動き検出を行っているが、拍動によるスフェアの動きとしては周期的な収縮・弛緩として現れるものもある。上記実施形態は画像の差分の大きさから動きを検出するので、このような収縮・弛緩であっても検出可能である。

また例えば、上記実施形態におけるスフェア抽出処理は、撮像された１つ１つの画像からスフェアのような比較的円形に近いオブジェクトを特定する処理であり、単独の静止画に対しても適用可能なものである。また、その対象物や利用目的も、心筋細胞スフェアやその拍動判定に限定されるものではなく、種々の細胞コロニーを対象として適用することが可能である。

また例えば、上記実施形態における拍動判定処理を実行するのに際して、予め画像内のスフェアを特定しておく必要があるが、スフェアを特定するための処理については上記に限定されるものではなく、他の方法で特定されたスフェアの情報を用いて、上記の拍動判定処理を実行してもよい。

またたとえば、上記実施形態では、撮像された各フレーム画像に対しスフェア抽出処理が実行され、フレームの間引きを行った上で拍動検出処理が実行される。これに代えて、例えば撮像された原画像に対しフレームの間引きを行い、間引き後のフレーム画像を用いてスフェア抽出処理および拍動検出処理が実行されてもよい。

以上、具体的な実施形態を例示して説明してきたように、この発明に係る画像処理方法では、同一の試料の同一範囲を互いに異なる撮像時刻に撮像した複数の原画像を取得し、複数の原画像各々について閉曲線の抽出、合成画像の作成およびオブジェクト分割を実行し、少なくとも一の合成画像に対するオブジェクト分割において、他の一の合成画像に対する距離変換の実行結果に基づき分水嶺法を実行するように構成されてもよい。このような構成によれば、経時的に動きのある細胞コロニーを画像間で相互に関連付けることが可能である。

この場合において例えば、複数の合成画像に対するオブジェクト分割において、単一の合成画像に対する距離変換の実行結果に基づき分水嶺法が実行されてもよい。このような構成によれば、複数の画像に含まれるオブジェクトを１つの画像中のオブジェクトを介して相互に関連付けることができる。

あるいは例えば、一の合成画像に対するオブジェクト分割において、対応する原画像の撮像時刻が当該合成画像と最も近い他の一の合成画像に対する距離変換の実行結果に基づき分水嶺法が実行されてもよい。このような構成によれば、撮像時間の差が小さく、したがってオブジェクトの変位が小さい画像間でオブジェクトの関連付けを行うことができるので、経時的なオブジェクトの移動を確実に捕捉することができる。

また、この発明に係る画像処理方法では、合成画像に対し距離変換を実行することで求められる距離の最大値が所定の閾値より小さいオブジェクトを除外して分水嶺法が実行されてもよい。このような構成によれば、面積の小さいオブジェクトや細いオブジェクト、外形において凹凸の多いオブジェクト等を除外することができる。

また、原画像の周縁部の全画素の濃度を背景に相当する濃度とみなしてもよい。このような構成によれば、原画像の端にかかるオブジェクトが閉じた輪郭を持たないものみなされて処理から除外されるのを回避することができる。

本発明は、細胞を撮像し画像を評価する生化学や医療の分野に適用可能であり、例えばｉＰＳ細胞やＥＳ細胞から分化培養される心筋細胞スフェア等の細胞コロニーを画像から検出する目的に特に好適である。

１ コンピュータ装置
２ 記録ディスク（記録媒体）
Ｉｂ ２値化画像
Ｉｆ フレーム画像（原画像）
Ｉｓ 合成画像
Ｓｐ スフェア（細胞コロニー）
Ｖｔｈ 閾値

Claims (8)

1. 培養された細胞コロニーを複数含む試料を撮像した原画像を取得する工程と、
   前記原画像から、所定のエッジ強度を有する閉曲線を抽出する工程と、
   抽出された前記閉曲線で囲まれた領域を示す画像と前記原画像を２値化した画像との論理積に相当する合成画像を作成する工程と、
   前記合成画像に含まれるオブジェクトに対し分水嶺法によるオブジェクト分割を実行する工程と
   を備える画像処理方法。
2. 同一の前記試料の同一範囲を互いに異なる撮像時刻に撮像した複数の前記原画像を取得し、前記複数の原画像各々について前記閉曲線の抽出、前記合成画像の作成および前記オブジェクト分割を実行し、
   少なくとも一の前記合成画像に対する前記オブジェクト分割において、他の一の前記合成画像に対する距離変換の実行結果に基づき前記分水嶺法を実行する請求項１に記載の画像処理方法。
3. 複数の前記合成画像に対する前記オブジェクト分割において、単一の前記合成画像に対する距離変換の実行結果に基づき前記分水嶺法を実行する請求項２に記載の画像処理方法。
4. 一の前記合成画像に対する前記オブジェクト分割において、対応する前記原画像の撮像時刻が当該合成画像と最も近い他の一の前記合成画像に対する距離変換の実行結果に基づき前記分水嶺法を実行する請求項２に記載の画像処理方法。
5. 前記合成画像に対し距離変換を実行することで求められる距離の最大値が所定の閾値より小さい前記オブジェクトを除外して前記分水嶺法を実行する請求項１ないし４のいずれかに記載の画像処理方法。
6. 前記原画像の周縁部の全画素の濃度を、背景に相当する濃度とみなす請求項１ないし５のいずれかに記載の画像処理方法。
7. コンピュータに、請求項１ないし６のいずれかに記載の画像処理方法の各工程を実行させるためのコンピュータプログラム。
8. 請求項７に記載のコンピュータプログラムを記録した、コンピュータ読み取り可能な記録媒体。

Images