(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-12-19
(45)【発行日】2022-12-27
(54)【発明の名称】画像処理方法、コンピュータプログラムおよび記録媒体
(51)【国際特許分類】
C12Q 1/04 20060101AFI20221220BHJP
G06T 7/00 20170101ALI20221220BHJP
G06T 7/254 20170101ALI20221220BHJP
C12M 1/00 20060101ALN20221220BHJP
C12M 1/34 20060101ALN20221220BHJP
【ＦＩ】
C12Q1/04
G06T7/00 630
G06T7/254 A
C12M1/00 A
C12M1/34 Z
【請求項の数】
11
(21)【出願番号】P 2018177298
(22)【出願日】2018-09-21
(65)【公開番号】P2020043842
(43)【公開日】2020-03-26
【審査請求日】2021-06-18
(73)【特許権者】
【識別番号】000207551
【氏名又は名称】株式会社ＳＣＲＥＥＮホールディングス
(74)【代理人】
【識別番号】100105935
【弁理士】
【氏名又は名称】振角 正一
(74)【代理人】
【識別番号】100136836
【弁理士】
【氏名又は名称】大西 一正
(72)【発明者】
【氏名】古田 智靖
【審査官】斉藤 貴子
(56)【参考文献】
【文献】特開２００８－０７６０８８（ＪＰ，Ａ）
【文献】特開２０１６－０６７２８７（ＪＰ，Ａ）
【文献】特開２０２０－０６４３４９（ＪＰ，Ａ）
【文献】国際公開第２０１５／０４１１７７（ＷＯ，Ａ１）
【文献】国際公開第２０１７／１５４３１８（ＷＯ，Ａ１）
【文献】特開２０１６－０４１０２３（ＪＰ，Ａ）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｃ１２Ｑ
Ｃ１２Ｍ
Ｇ０６Ｔ
ＣＡｐｌｕｓ／ＭＥＤＬＩＮＥ／ＥＭＢＡＳＥ／ＢＩＯＳＩＳ（ＳＴＮ）
ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）
(57)【特許請求の範囲】
【請求項1】
培養された心筋細胞のコロニーを含む試料を異なる撮像時刻に撮像した複数の原画像を取得する工程と、前記原画像のうち互いに異なる第１画像と第２画像とのそれぞれから前記コロニーに対応するコロニー領域を抽出する工程と、
前記第１画像と前記第２画像との間における前記コロニー領域の変位を検出しその検出結果に基づき前記コロニーの拍動の有無を判定する工程と
を備え、
前記コロニー領域のうち一の対象コロニー領域について、所定の探索範囲内における前記第１画像と前記第２画像との差分の大きさを前記対象コロニー領域の前記変位の大きさとして検出し、
前記探索範囲は、前記第１画像における前記対象コロニー領域の輪郭を所定幅に膨張させた環状の領域と、前記第２画像における前記対象コロニー領域の輪郭を前記所定幅に膨張させた環状の領域との論理和に相当する領域である、画像処理方法。
【請求項2】
培養された心筋細胞のコロニーを含む試料を異なる撮像時刻に撮像した複数の原画像を取得する工程と、前記原画像のうち互いに異なる第１画像と第２画像とのそれぞれから前記コロニーに対応するコロニー領域を抽出する工程と、
前記第１画像と前記第２画像との間における前記コロニー領域の変位を検出しその検出結果に基づき前記コロニーの拍動の有無を判定する工程と
を備え、
前記コロニー領域のうち一の対象コロニー領域について、所定の探索範囲内における前記第１画像と前記第２画像との差分の大きさを前記対象コロニー領域の前記変位の大きさとして検出し、
前記探索範囲は、前記第１画像における前記対象コロニー領域の内部領域と、前記第２画像における前記対象コロニー領域の内部領域との排他的論理和に相当する領域を少なくとも含む環状の領域である、画像処理方法。
【請求項3】
前記第１画像と前記第２画像との差分画像を２値化した２値化画像と、前記探索範囲に対応する探索範囲画像との論理積に相当する領域の面積に基づく値を、前記変位の大きさの指標値とする請求項１または２に記載の画像処理方法。
【請求項4】
前記第１画像における前記対象コロニー領域の輪郭長と前記第２画像における前記対象コロニー領域の輪郭長との平均値で前記面積を正規化した値を、前記面積に基づく値として前記変位の大きさの指標値とする請求項３に記載の画像処理方法。
【請求項5】
前記複数の原画像から選択した一の前記第１画像に対し、撮像時刻が互いに異なる複数の前記原画像を１つずつ順に前記第２画像に設定して、それぞれの組み合わせについて前記コロニー領域の前記変位量を検出する請求項１ないし４のいずれかに記載の画像処理方法。
【請求項6】
撮像時刻が互いに異なる複数の前記原画像を１つずつ順に前記第１画像に設定し、当該第１画像より後に撮像された前記原画像の一を前記第２画像として設定する請求項１ないし４のいずれかに記載の画像処理方法。
【請求項7】
前記コロニーの拍動の有無を判定する工程では、
前記第１画像と前記第２画像との間における前記対象コロニー領域の重心の移動方向を前記変位の方向とし、
前記対象コロニー領域について、所定の閾値より大きい変位が複数回検出され、かつそれらの変位に、方向が互いに反対であるものが含まれる場合に、当該対象コロニー領域に対応する前記コロニーについて拍動ありと判定する請求項５または６に記載の画像処理方法。
【請求項8】
２つの変位について、それぞれの変位方向を示す方向ベクトルのなす角が９０度より大きいときそれらの方向が反対であるとみなす請求項７に記載の画像処理方法。
【請求項9】
前記原画像から前記コロニー領域を抽出する工程は、前記原画像から所定のエッジ強度を有する閉曲線を抽出する工程と、
抽出された前記閉曲線で囲まれた領域を示す画像と前記原画像を２値化した画像との論理積に相当する合成画像を作成する工程と、
前記合成画像に含まれるオブジェクトに対し分水嶺法によるオブジェクト分割を実行する工程と
を含む請求項１ないし８のいずれかに記載の画像処理方法。
【請求項10】
コンピュータに、請求項１ないし９のいずれかに記載の画像処理方法の各工程を実行させるためのコンピュータプログラム。
【請求項11】
請求項１０に記載のコンピュータプログラムを記録した、コンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【0001】
この発明は、培養された心筋細胞のコロニーを撮像した画像から拍動の有無を判定する画像処理技術に関するものである。
【背景技術】
【0002】
再生医療や創薬などの分野においては、例えばｉＰＳ細胞（induced Pluripotent Stem Cell）やＥＳ細胞（Embryonic
Stem Cell）などの多能性幹細胞を分化培養し、心筋細胞を人工的に作製する技術が研究されている。良好に培養された心筋細胞では周期的な拍動運動がみられることから、培養状態を評価するために、心筋細胞のコロニーを含む試料を顕微鏡等で観察し拍動の有無を判定することが行われる。目視による判定は観察者の主観により結果が変動し、また多数のコロニーを観察することは観察者の負担となるため、試料を撮像した画像のデータ解析処理によって拍動の有無を自動的に検出する方法が提案されている。
Stem Cell）などの多能性幹細胞を分化培養し、心筋細胞を人工的に作製する技術が研究されている。良好に培養された心筋細胞では周期的な拍動運動がみられることから、培養状態を評価するために、心筋細胞のコロニーを含む試料を顕微鏡等で観察し拍動の有無を判定することが行われる。目視による判定は観察者の主観により結果が変動し、また多数のコロニーを観察することは観察者の負担となるため、試料を撮像した画像のデータ解析処理によって拍動の有無を自動的に検出する方法が提案されている。
【0003】
例えば特許文献１に記載の技術では、所定の時間間隔で撮像された時系列画像から動きベクトルが検出され、動きベクトルが検出された領域のうち内部構造の顕著な変化が見られるものが拍動領域として検出される。
【先行技術文献】
【特許文献】
【0004】
【文献】特開２０１４－０７５９９９号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
分化培養の初期段階にある試料では、幹細胞から分化した心筋細胞により形成される小さなコロニーが試料中に多数分散した状態にある。より具体的には、いくつかの心筋細胞が球状に集まってスフェアとも呼ばれる細胞コロニーを形成し、そのようなコロニーが試料内に散在する。それらのコロニーは互いに独立した周期で拍動し、また拍動しないものも含まれ得る。また拍動時のコロニーの動きは必ずしも収縮・弛緩ではなく、画像内では周期的な振動として現れることも多い。さらに、コロニーが高密度で存在する試料の画像では、コロニー同士が接触している場合が多い。
【0006】
上記従来技術は、比較的大きなコロニーにおける律動的な拍動を検出することに適している。しかしながら、上記したように小さなコロニーが多数分散し、個々のコロニーが個別に拍動する試料において、それを撮像した画像から個々のコロニーの拍動の有無を判定するのには必ずしも向いていないと考えられる。このため、このような小さなコロニーが多数ある場合に好適な拍動検出のための画像処理技術の確立が望まれる。
【0007】
この発明は上記課題に鑑みなされたものであり、培養された心筋細胞のコロニーを複数含む画像から、個々のコロニーの拍動の有無を判定するのに好適な画像処理技術を提供することを目的とする。
【課題を解決するための手段】
【0008】
この発明に係る画像処理方法は、上記目的を達成するため、培養された心筋細胞のコロニーを含む試料を異なる撮像時刻に撮像した複数の原画像を取得する工程と、前記原画像のうち互いに異なる第１画像と第２画像とのそれぞれから前記コロニーに対応するコロニー領域を抽出する工程と、前記第１画像と前記第２画像との間における前記コロニー領域の変位を検出しその検出結果に基づき前記コロニーの拍動の有無を判定する工程とを備えている。
【0009】
本発明の一の態様では、前記コロニー領域のうち一の対象コロニー領域について、所定の探索範囲内における前記第１画像と前記第２画像との差分の大きさを前記対象コロニー領域の前記変位の大きさとして検出し、前記探索範囲は、前記第１画像における前記対象コロニー領域の輪郭を所定幅に膨張させた環状の領域と、前記第２画像における前記対象コロニー領域の輪郭を前記所定幅に膨張させた環状の領域との論理和に相当する領域である。
【0010】
また、本発明の他の一の態様では、前記コロニー領域のうち一の対象コロニー領域について、所定の探索範囲内における前記第１画像と前記第２画像との差分の大きさを前記対象コロニー領域の前記変位の大きさとして検出し、前記探索範囲は、前記第１画像における前記対象コロニー領域の内部領域と、前記第２画像における前記対象コロニー領域の内部領域との排他的論理和に相当する領域を少なくとも含む環状の領域である。
【0011】
このように構成された発明では、異なる時刻に撮像された第１画像と第２画像とのそれぞれにおいてコロニー領域が抽出され、両画像間におけるコロニー領域の変位が検出される。この変位は、第１画像および第２画像において抽出されたコロニー領域の差分として検出される。ただし、コロニー領域同士の単純な差分ではなく、抽出されたコロニー領域の輪郭に基づき設定される探索範囲内での差分が求められる。
【0012】
詳しくは後述するが、上記した２つの態様のいずれにおいても、第１画像および第２画像でそれぞれ検出されたコロニー領域の輪郭の周囲に探索範囲が設定される。したがって、求められる差分の大きさは、コロニー領域の輪郭部分の動きの大きさを強く反映したものとなる。その一方、輪郭部分から遠く探索範囲に含まれないコロニー内部やコロニー外の変位は差分の値に影響を与えない。したがって、１つの探索範囲について求められる差分の大きさは、１つのコロニーの変位の大きさを表すものとなる。
【0013】
本願発明者の知見によれば、心筋細胞によるスフェアのような小さなコロニーについては、こうして求められる差分の大きさを当該コロニーの変位の大きさの指標とすることで試料中の各コロニーの変位を個別に、かつ的確に把握することができる。そして、こうして把握される変位の大きさから、当該コロニーの拍動の有無を判定することができる。原画像中に多数のコロニーがある場合であっても、各コロニーについてその輪郭周辺の変位のみが検出されるため、他のコロニーの挙動に影響されることなく拍動の有無を判定することが可能である。
【0014】
本発明に係る画像処理方法は、コンピュータをその実行主体とすることが可能なものである。この意味において、本発明は、コンピュータに上記処理を実行させるためのコンピュータプログラムとして実現することが可能である。また、当該コンピュータプログラムを記録した記録媒体として実現することも可能である。
【発明の効果】
【0015】
以上のように、本発明によれば、培養された心筋細胞のコロニーを複数含む画像から特定される各コロニーの輪郭の周辺における画像の差分が当該コロニーの変位の大きさの指標とされる。このため、画像に多数のコロニーが含まれる場合でも、他のコロニーの挙動に影響されることなく、個々のコロニーの拍動の有無を適切に判定することができる。
【図面の簡単な説明】
【0016】
【図1】本実施形態の画像処理を示すフローチャートである。
【図2】１つのフレーム画像を模式的に例示する図である。
【図3】前処理工程を示すフローチャートである。
【図4】スフェア抽出工程を示すフローチャートである。
【図5】画像端補正工程を説明するための図である。
【図6】エッジ検出処理後の画像の例を示す図である。
【図7】輪郭領域画像の一例を示す図である。
【図8】２値化画像の一例を示す図である。
【図9】論理積合成画像の一例を示す図である。
【図10】距離画像の一例を示す図である。
【図11】距離変換結果に基づく不要オブジェクト消去の概念を示す図である。
【図12】中心画像の一例を示す図である。
【図13】オブジェクト分割の概念を示す図である。
【図14】異なるフレーム間でのスフェアの対応関係を示す図である。
【図15】間引き処理の原理を示す図である。
【図16】動き検出工程を示すフローチャートである。
【図17】基準画像と比較画像との関係を示す図である。
【図18】基準画像と比較画像との差分を示す画像の例である。
【図19】探索範囲設定方法の第１の態様を示す図である。
【図20】探索範囲設定方法の第２の態様を示す図である。
【図21】収縮膨張画像と探索範囲画像との論理積合成画像を示す図である。
【図22】拍動判定工程を示すフローチャートである。
【図23】動き方向の判断例を示す図である。
【図24】この実施形態の画像処理を実行するコンピュータ装置の構成例である。
【発明を実施するための形態】
【0017】
＜概要の説明＞
以下、本発明に係る画像処理方法の一実施形態について説明する。この実施形態における画像処理は、培養された心筋細胞のコロニーを含む試料を撮像した画像から細胞コロニーを検出し、さらには各細胞コロニーの拍動の有無を判定するための処理である。例えば適宜の培地中でｉＰＳ細胞（induced Pluripotent Stem Cell）やＥＳ細胞（Embryonic
Stem Cell）などの多能性幹細胞から分化培養された心筋細胞のコロニーを含む試料を、本実施形態の処理対象物とすることができる。
以下、本発明に係る画像処理方法の一実施形態について説明する。この実施形態における画像処理は、培養された心筋細胞のコロニーを含む試料を撮像した画像から細胞コロニーを検出し、さらには各細胞コロニーの拍動の有無を判定するための処理である。例えば適宜の培地中でｉＰＳ細胞（induced Pluripotent Stem Cell）やＥＳ細胞（Embryonic
Stem Cell）などの多能性幹細胞から分化培養された心筋細胞のコロニーを含む試料を、本実施形態の処理対象物とすることができる。
【0018】
このような培養の初期段階では、幹細胞から分化した心筋細胞がいくつか集まって概ね球状のスフェアと呼ばれる小さなコロニーを形成する。試料には多数のスフェアが含まれ得るが、それらは互いに独立に拍動し、また中には拍動しないものもある。その他、未分化の細胞や細胞の代謝に伴い排出される老廃物等も試料には含まれ得る。以下に説明する本実施形態の画像処理は、このような画像から個々のスフェアを検出し、それらの拍動の有無を個別に判定することを目的とするものである。
【0019】
この画像処理は、心筋細胞スフェアを含む試料が、所定の時間間隔で繰り返し撮像する動画撮像機能を備えた撮像装置によって撮像された画像に対してなされる。撮像装置が撮像後の画像データに対する後処理の１つとして本画像処理を実行してもよく、また撮像装置あるいは適宜の記憶手段から撮像画像に対応する画像データを受け取ったコンピュータ装置が、本画像処理を実行してもよい。
【0020】
ここでは、汎用のコンピュータ装置が公知の撮像装置を制御して画像を取得し、該画像に対し本実施形態の画像処理を実行する態様を例示して説明する。撮像装置の構成については特に限定されず、培養された細胞等を含む試料を培地とともに動画撮像し、デジタル画像データとして出力する機能を有するものであればよい。例えば顕微鏡にデジタルビデオカメラを組み合わせた撮像装置を用いることができる。画像としては例えば明視野画像を好適に用いることができる。また、コンピュータ装置についても、例えばパーソナルコンピュータとして製品化されている一般的なハードウェア構成を有するものを利用可能である。
【0021】
図１は本実施形態の画像処理を示すフローチャートである。まず、図１を参照してこの実施形態における画像処理の各工程の概略について説明し、その後で各工程における処理の具体的内容について詳しく説明する。なお、この画像処理は、撮像した画像からスフェアを検出するスフェア検出処理（ステップＳ１０１～Ｓ１０３）と、検出された各スフェアの拍動を検出する拍動検出処理（ステップＳ１０５～Ｓ１０７）との２つの処理に大別することができる。また、これらの間に、それぞれの処理に適したフレームレートを実現するための間引き処理が設けられる。
【0022】
コンピュータ装置は、培地中で培養された心筋細胞のコロニーを含む試料を適宜のフレームレートで動画撮像する（ステップＳ１０１）。フレームレートとしては、心筋細胞の拍動周期として想定される周期よりは十分に短い撮像周期となることが好ましく、例えば一般的な動画撮像におけるフレームレートである３０ないし６０ＦＰＳ（Frame Per Second）程度とすることができる。また、撮像の継続時間としては、拍動する細胞コロニーの動きを１周期以上にわたり確実に捉えることのできる時間、例えば数秒とすることができる。
【0023】
図２は１つのフレーム画像を模式的に例示する図である。フレーム画像Ｉｆでは、略一様かつ低濃度（高輝度）である背景ＢＧの中に、心筋細胞が略球状に集まって形成されたスフェアＳｐに対応するオブジェクトが背景より高い濃度を有して多数分布している。この他に、未分化細胞、遊離細胞、細胞から排出された老廃物の塊等、スフェア以外の物体Ｄｅに対応するオブジェクトも画像内には存在する。なお、背景とオブジェクトとの濃度の関係が上記と逆である場合には、以下の説明における濃淡についても逆に読み替えればよい。なお、図を見やすくするために、全てのオブジェクトに符号を付すことをせず、代表的な一部のオブジェクトにのみ符号を付している。以下の各図においても同様である。
【0024】
図２に符号Ａ～Ｄを付して例示するように、複数のオブジェクトが画像内で互いに接し結合していることがある。このような結合が見られる理由としては、試料中で実際にオブジェクト同士が接触している場合と、試料内では離れているが画像内では接しているように見える場合とがある。
【0025】
こうして撮像された画像については、フレームごとに、画像の微小な乱れを是正するための前処理工程が実行される（ステップＳ１０２）。そして、前処理後のフレーム画像に対してスフェア抽出工程が実行される（ステップＳ１０３）。この処理は、画像内に点在するオブジェクトの中からスフェアとみられるものを抽出し、その輪郭位置を特定するための処理である。これにより、スフェアが画像中でどの領域を占めているかが特定される。画像に複数のスフェアが含まれる場合、それらの輪郭が個別に特定される。
【0026】
続いてフレーム間引きが行われる（ステップＳ１０４）。後述するように、本実施形態におけるスフェアの拍動検出は、フレーム画像間におけるスフェアの変位の大きさに基づいて行われる。ここで、試料の撮像はスフェアの拍動周期に対して十分に短い周期で実行されるが、撮像時刻の差が小さすぎる画像間ではスフェアの変位が小さくその検出が難しくなる。したがって全てのフレーム画像を用いてスフェアの動きを検出しようとすることは必ずしも現実的ではない。
【0027】
一方で、撮像時刻の差が大きすぎると周期的な変位を見落とすおそれがある。拍動周期を考慮した適宜の間隔で画像を間引くことで、スフェアの変位をより確実に捕捉することが可能となる。また、画像処理のための演算量を大きく低減させることにも寄与する。この目的のために、撮像されたフレーム画像の中から所定のサンプリング間隔で一部を選択するフレーム間引きが行われる。
【0028】
次に、画像間でのスフェアの動きが検出される（ステップＳ１０５）。上記したように、撮像時刻の異なる２つの画像間でのスフェアの変位から、当該スフェアの動きが検出される。そして、こうして検出された変位量から、各スフェアの拍動の有無が判定される（ステップＳ１０６）。画像内でのオブジェクトの変位は拍動によるものに限らない。検出された変位が拍動に起因するものであるか否かにより、拍動の有無が判定される。判定結果については適宜の形式で出力されユーザに提示される（ステップＳ１０７）。以上が本実施形態の画像処理の全体の流れである。以下、各処理の具体的内容につき順次説明する。
【0029】
＜スフェア検出処理＞
図３は前処理工程を示すフローチャートである。撮像された各フレームの画像に対し、後述の画像処理の条件を一定とするためのコントラスト補正（ステップＳ２０１）およびノイズ除去（ステップＳ２０２）が行われる。これらの処理としては公知のものを適用可能であり、ここでは説明を省略する。例えばノイズ除去処理としてはバイラテラルフィルタ等の空間フィルタリング技術を用いることができる。なお、撮像条件によってはこれらの前処理を省くことができる場合がある。
図３は前処理工程を示すフローチャートである。撮像された各フレームの画像に対し、後述の画像処理の条件を一定とするためのコントラスト補正（ステップＳ２０１）およびノイズ除去（ステップＳ２０２）が行われる。これらの処理としては公知のものを適用可能であり、ここでは説明を省略する。例えばノイズ除去処理としてはバイラテラルフィルタ等の空間フィルタリング技術を用いることができる。なお、撮像条件によってはこれらの前処理を省くことができる場合がある。
【0030】
次に、スフェア抽出工程（ステップＳ１０３）について説明する。スフェア抽出工程では、画像内に存在する各種オブジェクトの領域を特定し、それらの中からスフェアＳｐであるものを選択し、画像内で結合しているスフェアＳｐについては適切な境界で区分することで、個々のスフェアに対応する画像内のオブジェクトを特定する。この処理は撮像により得られたフレーム画像Ｉｆの全てについて実行される。
【0031】
図４はスフェア抽出工程を示すフローチャートである。最初に１つのフレーム画像を選択し（ステップＳ３０１）、選択されたフレーム画像Ｉｆに対し画像端補正を行う（ステップＳ３０２）。図２に示すように、画像内のオブジェクトのうちフレーム画像Ｉｆの端部付近にあるものでは、その一部が画像端にかかり、全体が画像内に収まっていない場合がある。後述する処理においては輪郭とみなせる閉曲線で囲まれた領域を抽出する工程があるが、このように画像端で輪郭が途切れているオブジェクトの抽出漏れを回避するために、画像端補正が実行される。
【0032】
図５は画像端補正工程を説明するための図である。図５（ａ）に示すように、フレーム画像Ｉｆの周縁部を取り囲むように、１画素分の幅を有する額縁状の領域Ｒｆを設ける。この領域の画素値（輝度値）は背景ＢＧと同程度、またはこれより低濃度（高輝度）となるように設定される。例えば濃度０（すなわち最大輝度の白色）とすることができる。この領域Ｒｆは、画像の端部に相当する画素の画素値を上記値に転換することによって設けられてもよく、また画像の端部の外側に上記画素値を有する仮想的な画素を付加することにより設けられてもよい。
【0033】
処理後のフレーム画像Ｉｆでは、図５（ｂ）に示すように、画像端部の画素は全て背景濃度と同じまたはこれより低濃度となっている。このため、画像端にかかる（例えば図中に符号Ａ、Ｂで示す）オブジェクトも、画像内で閉じた輪郭を有することになる。これにより、輪郭の途切れに起因するオブジェクトの抽出漏れが回避される。なお図５（ｂ）では、画像端の様子を見やすくするためにフレーム画像Ｉｆの周縁部を点線で示している。
【0034】
図４に戻って、こうして画像端が補正されたフレーム画像Ｉｆに対して、輝度値の変化が大きな箇所を検出するエッジ検出処理が実行される（ステップＳ３０３）。エッジ検出処理としては公知の方法、例えばSobelフィルタに代表される空間フィルタ処理を適用することができる。
【0035】
図６はエッジ検出処理後の画像の例を示す図である。図２に示すフレーム画像Ｉｆに対応するエッジ検出処理後の画像（エッジ画像）Ｉｅでは、画像内のオブジェクトと背景との境界、つまりオブジェクトの外周部に相当する部分が強調されエッジとして検出される。図に符号Ａを付して例示するように、オブジェクトの内部テクスチャに起因するエッジも含まれ得る。また、図に符号Ｂを付して例示するように、予め画像端補正を行っているため、画像端にかかるオブジェクトについても閉じた輪郭が形成される。なお、検出された輪郭の微小な欠落を補償するために適宜の補正処理、例えば膨張収縮処理がさらに実行されてもよい。
【0036】
こうして検出されたエッジが閉曲線をなすとき、該閉曲線は画像中のオブジェクトの輪郭である可能性が高い。そこで、エッジ画像Ｉｅにおいて閉曲線で囲まれた輪郭の内部領域を抽出しこれを塗りつぶすことで（ステップＳ３０４）、輪郭の内部と外部とを区分する。ここでは、塗りつぶされた輪郭内部の領域を「輪郭領域」と称し、塗りつぶし後の画像を「輪郭領域画像」と称することとする。
【0037】
図７は輪郭領域画像の一例を示す図である。この輪郭領域画像Ｉｒは、図６のエッジ画像Ｉｅにおいて閉曲線で囲まれた領域を淡色（白）で、それ以外の領域を濃色（黒）で塗り分けた画像に相当する。なお、ここでは一般的な画像処理の例に倣いこのような塗り分けとしているが、濃淡が逆であってもよい。以下の処理においても同様である。
【0038】
輪郭領域画像Ｉｒでは、閉曲線で囲まれたオブジェクトの内部が全て塗りつぶされるため、オブジェクトの内部テクスチャに対応する情報は失われる。また、図２に示す元のフレーム画像Ｉｆと対比すると、図７に符号Ａ、Ｂを付して例示するように、本来は背景部分であるが、周囲全てをオブジェクトによって取り囲まれているために塗りつぶされてしまう領域も生じる。すなわち、本来はオブジェクトの外部とされるべき背景領域がオブジェクトの内部として扱われることとなってしまう。
【0039】
この問題を回避するための処理が次に実行される。すなわち、元のフレーム画像Ｉｆを適宜の閾値により２値化した画像が生成され（ステップＳ３０５）、さらにこの２値化画像と輪郭領域画像Ｉｒとを論理積合成した合成画像が生成される（ステップＳ３０６）。２値化されるフレーム画像Ｉｆとしては、前処理後のものが用いられることが好ましい。なお２値化のための閾値については固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としてもよい。
【0040】
図８は２値化画像の一例を示す図である。より具体的には、図８は、図２のフレーム画像Ｉｆに前処理を施した画像を適宜の閾値で２値化し、閾値より高濃度の領域を淡色（白）で、より低濃度の領域を濃色（黒）で塗り分けることで得られる２値化画像Ｉｂを示している。２値化画像Ｉｂでは、オブジェクトの輪郭部分およびその内部の高濃度領域が白色で表される一方、元のフレーム画像Ｉｆにおいて低濃度であった、背景領域やオブジェクト内の一部の低濃度領域が黒色で表される。図８に符号Ａ、Ｂを付して例示するように、周囲をオブジェクトで囲まれた背景部分も黒色で表される。なお、画像Ｉｂの左上隅および右上隅の本来は背景である部分に現れている白色の領域は、背景濃度のムラが２値化の結果に影響を及ぼす可能性があることを象徴的に示したものである。
【0041】
図９は論理積合成画像の一例を示す図である。図７に示す輪郭領域画像Ｉｒと図８に示す２値化画像Ｉｂとを論理積合成すると、図９に示す合成画像Ｉｓのように、元の両画像においてともに白色であった領域は白色で、それ以外の領域は黒色で表されることになる。これにより、図７の輪郭領域画像Ｉｒでは白色に塗りつぶされていた、オブジェクトで囲まれた背景領域が、他の背景部分と同じ黒色で表されることになる。一方で、２値化画像Ｉｂに含まれていたオブジェクトの内部テクスチャに起因する黒色領域は合成画像Ｉｓにおいても保存されている。
【0042】
したがって、図９に示す合成画像Ｉｓは、オブジェクトに囲まれた背景を正しく背景として認識するという意味において、またオブジェクトの内部テクスチャに関する情報を保存するという意味において、図７の輪郭領域画像Ｉｒよりも有効なものであるといえる。また、合成画像Ｉｓは、明確なエッジとして検出された輪郭によって、背景の濃度ムラの影響を受けずオブジェクトと背景との境界が画定されるという意味において、図８の２値化画像Ｉｂよりも有効なものであるといえる。つまり、合成画像Ｉｓは、画像中のオブジェクトと背景とを区分するという目的において、輪郭抽出に基づく方法の持つ利点と、２値化に基づく方法の持つ利点とを併せ持ったものである。
【0043】
ここまでの処理により、画像内でオブジェクトが占める領域と背景領域とが的確に区別される。しかしながら、前述したように、画像内のオブジェクトには、本実施形態の処理対象であるスフェア以外の不要なオブジェクトも含まれ、また複数のオブジェクトが画像内で接触、結合している場合がある。このため、不要オブジェクトの除去および結合したオブジェクトの分割のための処理が必要である。
【0044】
この目的のために、画像内の各オブジェクトの構造解析を行う。具体的には、ステップＳ３０６で生成された合成画像Ｉｓに対し距離変換が実行される（ステップＳ３０７）。２次元２値画像における距離変換は、汎用の画像処理ソフトウェアでも実行可能な一般的な画像処理であり、画像内の各画素を当該画素から最も近い黒色画素までの距離によって表すというものである。
【0045】
この実施形態において処理対象となる合成画像Ｉｓでは、背景に対応する画素において距離値０、背景画素に隣接しオブジェクトの輪郭に対応する画素において距離値１を取り、オブジェクト内部の画素ほど背景領域から遠ざかるため距離値は大きくなる。画像内の各画素の画素値を距離値に置き換え、例えば距離値が大きい画素ほど高い輝度値によって表すようにすれば、背景からの距離を濃淡で表した距離画像が得られる。
【0046】
図１０は距離画像の一例を示す図であり、この距離画像Ｉｄは図９に示す合成画像Ｉｓを距離変換した結果を示すものである。距離変換により、合成画像Ｉｓは、背景からの距離が各画素の輝度値に反映された距離画像Ｉｄが得られる。距離画像Ｉｄは、背景部分が暗く、オブジェクトの周縁部から内部に向かうにつれて高輝度となるグレースケール画像である。したがって、面積が大きく円形に近いオブジェクトでは中央部の距離値が大きく、距離画像では輝度が高く（明るく）なる。一方、面積の小さいオブジェクトや細いオブジェクト、周囲の凹凸の大きいオブジェクト等では中央部の距離値は小さく、したがって輝度値が低く（暗く）なる。
【0047】
処理対象物であるスフェアＳｐは、画像内で比較的大きく円形に近い形状を有する。したがって当該オブジェクトの内部には比較的大きな距離値を有する画素が存在する。一方、浮遊細胞等の小さなオブジェクトでは、オブジェクト内の各画素の距離値はより小さい。このことを利用して、画像内の各オブジェクトはスフェアＳｐとそれ以外の不要オブジェクトとに区別され、不要オブジェクトは距離画像Ｉｄおよび合成画像Ｉｓから消去される（ステップＳ３０８）。
【0048】
図１１は距離変換結果に基づく不要オブジェクト消去の概念を示す図である。図１１（ａ）に示すように、距離画像Ｉｄの各画素の位置を座標（Ｘ，Ｙ）で表し、各画素の距離値を座標平面に対しプロットした場合を考える。なお、図ではＸ座標およびＹ座標を１つの座標軸で表しているが、より厳密にはＸＹ座標平面に対し距離値の軸が直交する形態となる。
【0049】
このように表現したとき、距離値０の領域は背景に対応し、背景内で孤立した１つのオブジェクトは、０より大きい距離値を有する画素が連続する領域として表される。図１１（ａ）の例では、２つのオブジェクトＯｂｊ１、Ｏｂｊ２が互いに孤立して存在する。比較的サイズが大きく、また凹凸の少ないまとまった形状を有するオブジェクトＯｂｊ１では距離値のピークが大きく、より小さい細いオブジェクトＯｂｊ２ではピークが小さくなる。このことから、距離値に対し適宜の閾値Ｖｔｈを設定し、距離値のピーク値が閾値Ｖｔｈより大きいオブジェクトＯｂｊ１についてはスフェアとみなして画像に残す一方、距離値のピーク値が閾値Ｖｔｈより小さいオブジェクトＯｂｊ２については不要オブジェクトとみなして画像から消去する。
【0050】
このようにすれば、比較的大きなオブジェクト、円形に近いオブジェクト等は画像に残る一方、小さなオブジェクトや細いオブジェクト、周囲の凹凸の激しいオブジェクト等については画像から消去される。図１１（ｂ）は不要オブジェクト消去後の距離画像Ｉｄａの一例を示している。符号Ａ～Ｄを付した位置では、図１０に現れていた小さなオブジェクトや細長いオブジェクトが現れていない。すなわち、このようなオブジェクトは消去されている。画像の状態に応じた的確な不要オブジェクト除去を可能とするために、上記閾値Ｖｔｈについては固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としておくことが好ましい。
【0051】
このようにして、不要オブジェクトを除外した距離画像とこれに対応する中心画像Ｉｃａとを作成することができる（ステップＳ３０９）。ここでいう中心画像Ｉｃａは、距離画像Ｉｄａにおいて、距離値が極大となり、かつその極大値が所定値より大きい画素の位置を表す画像である。
【0052】
図１２は中心画像の一例を示す図である。より具体的には、図１２（ａ）は図１１（ｂ）に示す距離画像Ｉｄａに対応する中心画像Ｉｃａである。また、図１２（ｂ）は、図１２（ａ）の中心画像Ｉｃａを、図９に示す合成画像Ｉｓに重畳した参考図である。なお、図１２（ｂ）においては、中心画像Ｉｃａにより示される各点の視認性を高めるため、その位置を黒丸印で示す。
【0053】
これらの図からわかるように、また図１２（ｂ）に符号Ａ，Ｂを付して例示したように、円形に近いオブジェクトにあっては、このような画素は当該オブジェクトの中心に近い位置にある。この意味において、中心画像は各オブジェクトの中心位置を表しているといえる。そこで、中心画像Ｉｃａに現れる点を「中心点」と称することとする。図１２（ｂ）の画像には中心点を含まないオブジェクトも含まれているが、これらは不要オブジェクトとみなされるものであり、中心点の導出対象から除外されている。
【0054】
図１２（ｂ）に符号Ｃ，Ｄを付して例示するように、複数のオブジェクトが結合しているとみられる複雑な形状のオブジェクトや、符号Ｅを付した顕著なくびれを有するオブジェクトでは、１つのオブジェクトに複数の中心点が存在する。言い換えれば、このように複数の中心点が存在するオブジェクトは、中心点の数に対応する複数のオブジェクトが結合したものと考えられる。
【0055】
そこで、複数の中心点が存在するオブジェクトについては、それらの中心点の間に適宜の境界線を設定することで、それぞれ中心点を１つずつ有するオブジェクトに分割する。こうすることで、結合したオブジェクトを分割することができる。このような考え方は、「分水嶺（Watershed）法」と呼ばれるオブジェクト分割のための画像処理のものに他ならない。すなわち、距離変換と不要オブジェクト消去とにより得られた距離画像Ｉｄａおよび中心画像Ｉｃａの情報に基づき、合成画像Ｉｓに対し分水嶺法によるオブジェクト分割を実行することで（ステップＳ３１０）、画像内で結合するオブジェクトが分割される。
【0056】
図１３はオブジェクト分割の概念を示す図である。図１３（ａ）に模式的に示すように、２つの中心点Ｃ１、Ｃ２を有するオブジェクトＯｂｊ３は、それらの中心点Ｃ１、Ｃ２の間に、破線で例示する境界線を設定することで２つに分割される。この境界線を適切に設定するための画像処理アルゴリズムの１つが分水嶺法である。３以上の中心点を含むオブジェクトについても同様にすることができる。
【0057】
図１３（ｂ）はオブジェクト分割後の画像の例であり、より具体的には、図１１（ｂ）の距離画像Ｉｄａおよび図１２（ａ）の中心画像Ｉｃａの情報に基づく分水嶺法により、図９に示す合成画像Ｉｓをオブジェクト分割して得られたオブジェクトの輪郭を示す画像Ｉｏである。ただし、予め合成画像Ｉｓから不要オブジェクトは消去されているものとする。符号Ａ～Ｃを付して例示するように、合成画像Ｉｓでは結合していた複数のオブジェクトが分離され、また符号Ｄを付して例示するように、顕著なくびれを有するオブジェクトもそのくびれ部分から２つに分割される。こうして得られた１つ１つのオブジェクトは、それぞれ１つのスフェアとみなすことができる。
【0058】
このようにして１つのフレーム画像Ｉｆから抽出されたスフェアＳｐについてはインデックスが作成される（ステップＳ３１１）。例えば各スフェアＳｐに付した固有の識別記号と、その代表点（例えば中心）の位置を表す情報とを関連付けたインデックスを作成することができる。その他にも、例えば各スフェアＳｐの面積や輪郭の形状、周囲長などを表す情報をインデックスに含めてもよい。
【0059】
こうして１つのフレーム画像Ｉｆについてスフェア抽出が終了する。次に処理すべきフレーム画像がある場合には（ステップＳ３１２においてＹＥＳ）、ステップＳ３０１からの処理を繰り返すことで、各フレーム画像が順次処理される。最後のフレーム画像まで処理が終了すれば（ステップＳ３１２においてＮＯ）、フレーム抽出処理は完了する。
【0060】
上記したステップＳ３０１～Ｓ３１１の処理を各フレーム画像Ｉｆについて一律に実行することで、当該フレーム画像Ｉｆ中のスフェアＳｐを抽出することが可能である。この場合、複数のフレーム画像からそれぞれ抽出されるスフェアの間に何ら関連付けはなされない。しかしながら、実際には、小さな時間差で撮像された２つの画像でそれぞれ抽出されるスフェアは、ほぼ１対１の関係で関連付けることができる。言い換えれば、１つの画像に含まれるスフェアの各々は、多少の位置や形状の変化はあるとしても、他の１つの画像にもほぼ同様の位置に同じような形状で存在しているはずである。
【0061】
スフェアの拍動を検出するという目的に鑑みれば、１つのスフェアの位置が経時的にどう変化するかを追跡する必要がある。したがって、１つの画像に含まれるスフェアが他の画像中のどのスフェアに対応するかを把握しなければならない。このような観点から、第２フレーム以降の画像に対する処理は以下のように改変することができる。
【0062】
すなわち、改変されたスフェア抽出処理において、第２フレームの画像に対する処理では、ステップＳ３０９として記載された距離画像および中心画像の生成が省略される。そして、ステップＳ３１０の分水嶺法によるオブジェクト分割においては、第１フレームの画像に対する処理で生成された距離画像および中心画像の情報が使用される。つまり、第２フレームの画像中の各オブジェクトについては、その中心点が第１フレーム画像中のオブジェクトと同じ位置にあるものとみなされる。このようにする理由につき以下に説明する。
【0063】
図１４は異なるフレーム間でのスフェアの対応関係を示す図である。図１４（ａ）は画像中で孤立したスフェアの例を示している。なお、図１４（ａ）はフレーム画像の一部のみを示しており、また点線による格子は画像内でのスフェアの位置を把握しやすくするために加えた仮想的なものである。図左に示す第Ｎフレームの画像に含まれるスフェアＳｐは、図右に示す第（Ｎ＋１）フレームの画像では少し変位した位置に現れる。このため、第（Ｎ＋１）フレーム画像におけるスフェアＳｐの中心点Ｃ（Ｎ＋１）は、第Ｎフレーム画像におけるスフェアＳｐの中心点Ｃ（Ｎ）の位置から少し移動していることがある。
【0064】
しかしながら、短い撮像間隔の間での変位量は限定的であり、第Ｎフレーム画像における中心点Ｃ（Ｎ）は、第（Ｎ＋１）フレーム画像におけるスフェアＳｐの内部に含まれている、またはそれに極めて近い位置にあると予想される。したがって、第（Ｎ＋１）フレーム画像において中心点Ｃ（Ｎ）を含むスフェアが、または第（Ｎ＋１）フレーム画像において中心点Ｃ（Ｎ）に最も近い位置にあるスフェアが、第Ｎフレーム画像におけるスフェアと同一であるといえる。このことから、異なるフレーム間で対応するスフェアを相互に関連付けることができる。
【0065】
図１４（ｂ）は画像中で結合したスフェアの例を示している。第Ｎフレームの画像に２つのスフェアＳｐ１，Ｓｐ２が含まれ、これらが画像内で結合しているとする。スフェアＳｐ１，Ｓｐ２それぞれの中心点をＣ１（Ｎ），Ｃ２（Ｎ）とする。また、このときにオブジェクト分割により画定されるスフェア境界線を破線Ａによって表す。
【0066】
第（Ｎ＋１）フレームの画像において、図に点線で示すように、スフェアＳｐ２が右方向に変位したとする。このときのスフェアＳｐ２の中心点をＣ２（Ｎ＋１）とすると、スフェアＳｐ２の動きに対応して中心点Ｃ２（Ｎ＋１）も右方向に変位していると考えられる。一方、スフェアＳｐ１はほとんど動かず、第（Ｎ＋１）フレームにおける中心点Ｃ１（Ｎ＋１）は、第Ｎフレームにおける中心点Ｃ１（Ｎ）とほぼ同じである。
【0067】
第（Ｎ＋１）フレームの画像において、検出される中心点Ｃ１（Ｎ＋１），Ｃ２（Ｎ＋１）の位置に基づきオブジェクト分割を行ったとき、破線Ｂで示すように、中心点間の距離の増大に起因して境界線も右方向にシフトすることになる。これにより、スフェアＳｐ１が右方向に広がったように見える。しかしながら、実際にはスフェアＳｐ２が動いたのであって、スフェアＳｐ１が占める領域はほぼ変わらない。このようなオブジェクト分割の結果と実態との乖離を抑えるには、第（Ｎ＋１）フレームにおいても、第Ｎフレームで特定された中心点Ｃ１（Ｎ），Ｃ２（Ｎ）の情報を用いるようにすればよい。
【0068】
第（Ｎ＋２）フレーム以降の画像についてオブジェクト分割を行う際に、どのフレームで特定された中心画像の情報を用いるかについてはいくつかの考え方があり得る。例えば、全てのフレームの処理に対して単一のフレーム、例えば第１フレームにおける中心画像の情報を適用することが考えられる。また、第（Ｎ＋１）フレームの処理に対してその直前に撮像された第Ｎフレームにおける中心画像の情報を適用することが考えられる。さらに、一定のフレーム数ごとに中心画像を更新し、次の更新までは同じ情報を用いて各フレームの処理を行うようにしてもよい。
【0069】
このように、一のフレーム画像に対するオブジェクト分割に際し、異なる時刻に撮像されたフレーム画像の距離情報および中心情報を用いて分水嶺法を実行することにより、それらのフレーム画像間において対応するスフェアの関連付けを行うことができる。例えば各フレーム画像で検出されたスフェアにつきインデックスを作成する際に、対応するスフェアに共通の識別記号を付すようにすれば、画像間でのスフェアの対応関係を明確に記録することができる。
【0070】
以上のように、本実施形態のスフェア検出処理では、連続撮像された各フレーム画像からスフェアを検出するのに際して、次のような要素技術を含んでいる。まず、画像端にかかるオブジェクトが除外されるのを防止するために、画像の周縁部を背景とみなせる濃度に設定し、画像端にかかるオブジェクトが閉じた輪郭で囲まれるようにする。
【0071】
また、エッジ検出により元のフレーム画像から抽出されたオブジェクトの輪郭に基づく輪郭領域画像と、フレーム画像を濃度（輝度）に基づいて２値化した２値化画像との論理積合成画像が、画像においてスフェア等のオブジェクトと背景との境界を示すものとされる。これにより、オブジェクトに取り囲まれた背景領域が正しく背景として認識されるようになる。
【0072】
そして、距離変換の結果に基づきスフェアでない可能性の高いオブジェクトを消去し、残るオブジェクトに対し分水嶺法によるオブジェクト分割を実行することで、画像内で結合したスフェアをそれぞれ個別のオブジェクトとして特定することができる。このとき、撮像時刻の異なる画像における距離変換の結果を用いてオブジェクト分割を行うことで、経時的に動きのあるスフェアを画像間で関連付けることができ、また動きのないスフェアと動きのあるスフェアとが結合している場合でも、それらを的確にかつ安定的に分割することができる。
【0073】
＜間引き処理＞
図１５は間引き処理の原理を示す図である。図１５（ａ）に示すように、試料を所定のフレームレートで連続撮像することにより、一定の時間間隔で撮像された複数の原画像Ｉｇ１，Ｉｇ２，…が取得される。例えばフレームレートが６０ＦＰＳであるとき、時系列において隣り合う、つまり撮像時間間隔が最も小さい２つの原画像間の撮像時間間隔は（１／６０）秒である。一方で、心筋細胞スフェアの拍動周期は概ね１秒程度であるため、隣り合う２つの原画像の間では拍動に起因するスフェアの変位は極めて小さい。拍動に起因するスフェアの変位を効率よくかつ確実に検出するためには、拍動周期を考慮して画像の間引きを行うことが有効である。
図１５は間引き処理の原理を示す図である。図１５（ａ）に示すように、試料を所定のフレームレートで連続撮像することにより、一定の時間間隔で撮像された複数の原画像Ｉｇ１，Ｉｇ２，…が取得される。例えばフレームレートが６０ＦＰＳであるとき、時系列において隣り合う、つまり撮像時間間隔が最も小さい２つの原画像間の撮像時間間隔は（１／６０）秒である。一方で、心筋細胞スフェアの拍動周期は概ね１秒程度であるため、隣り合う２つの原画像の間では拍動に起因するスフェアの変位は極めて小さい。拍動に起因するスフェアの変位を効率よくかつ確実に検出するためには、拍動周期を考慮して画像の間引きを行うことが有効である。
【0074】
例えば１秒の周期で振動するスフェアの変位を検出するためには、フレーム間の撮像時間間隔は１／６秒ないし１／３秒程度が好ましい。この場合、原画像を１／１０ないし１／２０程度に間引くことになる。この間引きの度合いを以下では「間引き率」と称する。種々の試料に柔軟に対応するために、間引き率については固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としておくことが好ましい。なお、この場合の原画像Ｉｇ１，Ｉｇ２，…としては、例えばスフェア検出工程における前処理後のフレーム画像を用いることができる。
【0075】
図１５（ｂ）に示すように、原画像Ｉｇ１，Ｉｇ２，…から、間引き率に応じた一定のサンプリング間隔でサンプリング抽出した画像を時系列に沿って再配置することで、間引きが実現される。間引き後の画像群Ｉｆ１，Ｉｆ２，Ｉｆ３，…は、元のフレームレートに間引き率の逆数を乗じた値を実効的なフレームレートとする連続画像とみなすことができる。以下の処理は、これらの各画像を実質的なフレーム画像として実行される。そこで、以下の処理の説明において単に「フレームレート」および「フレーム画像」というとき、これらは間引き後のものを指すこととする。
【0076】
＜拍動検出処理＞
次に、上記のようにして検出された各スフェアについての拍動検出処理について説明する。本実施形態の拍動検出処理は、図１に示すように、動き検出工程（ステップＳ１０５）および拍動判定工程（ステップＳ１０６）を備えている。これらの工程は、画像中のスフェアの各々について実行される。
次に、上記のようにして検出された各スフェアについての拍動検出処理について説明する。本実施形態の拍動検出処理は、図１に示すように、動き検出工程（ステップＳ１０５）および拍動判定工程（ステップＳ１０６）を備えている。これらの工程は、画像中のスフェアの各々について実行される。
【0077】
図１６は動き検出工程を示すフローチャートである。また、図１７は基準画像と比較画像との関係を示す図である。この動き検出工程においては、間引き後のフレーム画像Ｉｆ１，Ｉｆ２，…から選択した１つの画像を基準画像とし（ステップＳ４０１）、これと撮像時刻が近いが異なる他の１つの画像を比較画像として（ステップＳ４０２）、両画像間におけるスフェアの変位を求め、その結果に基づき両画像の撮像時刻の間でのスフェアの動きを検出する。まず、基準画像と比較画像との関係につき、図１７を参照して説明する。
【0078】
スフェアの拍動周期に応じて適宜に間引かれたフレーム画像Ｉｆ１，Ｉｆ２，…における動き検出では、まず最初の画像Ｉｆ１が基準画像とされる。そして、これより後に撮像された画像Ｉｆ２，Ｉｆ３，Ｉｆ４を順次比較画像に設定し、それぞれ基準画像Ｉｆ１との間でスフェアの変位が求められる。
【0079】
時系列において隣り合う画像間での比較により、フレーム画像間の撮像時間間隔に等しい時間におけるスフェアの動きがわかる。ただし、スフェアの拍動周期は必ずしも一定ではなく、またスフェアごとに独立である。このため、基準画像と比較画像との撮像時間の差が、拍動に起因するスフェアの動きを捕捉するのに最適なものでない場合があり得る。例えば、拍動するスフェアが２つの画像の間で偶然に同じ位置にあるということが起こり得る。このことに起因する検出漏れを回避するために、比較画像を順次入れ替え、数種類の時間間隔でスフェアの変位を検出する。そして、基準画像と比較画像とのいずれかの組み合わせにおいてスフェアの変位が認められると動きありと判断するようにすれば、より確実な動き検出が可能となる。
【0080】
ここでは１つの基準画像に対し３つの比較画像を設定することとしているが、比較画像として用いる画像の数はこれに限定されず任意である。ただし、比較画像の数を増やすと演算処理量も増大し実行装置の負荷が大きくなる、そのため、試料の実状に応じて比較画像の数と間引き率とを適切に設定することで、上記のような検出漏れの防止を図ることが望まれる。
【0081】
基準画像Ｉｆ１に対し３つの比較画像を適用した後、次に基準画像を画像Ｉｆ２に切り替え、これに対しても同様に３つの画像Ｉｆ３，Ｉｆ４，Ｉｆ５を順次比較画像に設定してスフェアの変位が求められる。その後も同様に、基準画像を画像Ｉｆ３，Ｉｆ４，…と順番に切り替えながら、その都度３つの比較画像を順次設定して変位を求める。図１６のフローチャートでは、ループ処理によってこれが実現されている。このようにすることで、種々の周期で拍動するスフェアの動きを漏れなく検出することが可能となる。
【0082】
図１６に戻り動き検出工程の処理内容についての説明を続ける。一の基準画像および一の比較画像が選出されると（ステップＳ４０１，Ｓ４０２）、両画像間のコントラストの違いに起因する検出誤差を低くするために、コントラスト補正が実行される（ステップＳ４０３）。補正処理としては公知のものを適用可能である。また、補正が必要なければ省くことも可能である。
【0083】
続いて、基準画像と比較画像との差分画像が生成される（ステップＳ４０４）。差分画像はさらに２値化される（ステップＳ４０５）。両画像間でスフェアに動きがなければ画像内容に変化はなく、差分はゼロとなるはずである。一方、変位が大きければ差分も大きくなる。したがって、差分画像を２値化した画像は変位の大きさを定量的に表すための情報を有している。
【0084】
図１８は基準画像と比較画像との差分を示す画像の例である。互いに対応するスフェアＳｐを含む基準画像Ｉ１と比較画像Ｉ２とを考える。図において点線による格子は画像内でのスフェアの位置を把握しやすくするために加えた仮想的なものである。基準画像Ｉ１と比べると、比較画像Ｉ２におけるスフェアＳｐは、その外形形状を変化させながらわずかに右方向に変位している。これらの差分を求めると画像Ｉ３が得られる。差分画像Ｉ３では、スフェアＳｐの外周部の変位に起因する概ね環状の像と、その内部に含まれるスフェア内部のテクスチャの変化に起因する像とが現れている。実際の差分画像は多階調画像であり、定量的な取り扱いを容易にするために２値化することで、２値化画像Ｉ４が得られる。
【0085】
２値化画像Ｉ４に対し、公知の収縮膨張処理が実行される（ステップＳ４０６）。これにより、２値化画像Ｉ４から微細な斑点状のノイズ成分が除去される。なお、以下ではノイズ除去後の２値化画像を「収縮膨張画像」と称し符号Ｉ４ａにより表すこととする。この２値化画像Ｉ４ａは、スフェアの外周部の変位に起因する差分像とスフェア内部のテクスチャ変化に起因する差分像とを含む。このうちテクスチャ変化に起因する像は、スフェアの変位量を定量的に求めるに当たっては誤差要因となるものである。その理由は以下の通りである。
【0086】
内部テクスチャの変化は、変位に直接関係のない細胞の活動によっても生じる。このため、スフェアが変位せず内部テクスチャのみが変化する場合がある。このような変化が変位と誤認されることがないようにする必要がある。また、仮に内部テクスチャに変化がないとしても、スフェア全体が変位することによって、同一位置の画素同士の差分を求めれば有意な差が現れることとなる。この差の大きさはスフェアの内部構造に依存するから、差分の値は変位量のみならずスフェア内部構造にも影響を受けていることになる。より具体的には、同量の変位であってもスフェアの内部構造によって差分の大きさが変わる。
【0087】
したがって、スフェアの変位のみを定量的に求めるためには、スフェアの内部構造の影響を排除し、スフェア外周部の変位に起因する差分のみを抽出する必要がある。このために、本実施形態では、スフェア外周部の近傍に探索範囲を設定し（ステップＳ４０７）、この探索範囲内での画像間の差分を有意な値として求めるという手法が採用されている。以下、探索範囲設定方法の２つの態様につき説明するが、いずれの方法が適用されてもよい。
【0088】
図１９は探索範囲設定方法の第１の態様を示す図である。より具体的には、図１９（ａ）は探索範囲設定方法の第１の態様を示すフローチャートであり、図１９（ｂ）は処理過程における画像の例を模式的に示す図である。図１９（ａ）に示すように、基準画像中のスフェアおよび比較画像中のスフェアのそれぞれにつき輪郭描画が行われる（ステップＳ５０１、Ｓ５０２）。各フレーム画像におけるスフェアの輪郭は、先のスフェア抽出処理で既に求められている。こうして求められた輪郭情報に基づき輪郭描画が行われるが、その線幅は１画素より大きな値に設定される。つまり、輪郭を膨張させた状態で描画する。
【0089】
輪郭描画の結果、図１９（ｂ）に示すように、基準画像および比較画像それぞれにおいて特定されたスフェアの輪郭（点線）を所定の線幅に膨張させて描画した輪郭画像Ｉ５、Ｉ６が得られる。これらを論理和合成すると（ステップＳ５０３）、図１９（ｂ）下部に示すように、合成された画像Ｉ７には、基準画像および比較画像それぞれにおけるスフェアの輪郭を共に含み、かつそれらの輪郭の周囲まで広がる環状の帯状領域が現れる。この帯状領域が探索範囲Ｒｓとされる。ここでは画像Ｉ７を第１の態様における「探索範囲画像」と称することとする。
【0090】
外周部の位置の変化を確実に差分として捕捉するためには、膨張させる輪郭の線幅は太い方が望ましい。その一方で、そうして探索範囲を広くすると、スフェアの外周部の変化以外の変化、例えば外周部に近い領域の内部テクスチャの変化も検出されてしまうことになる。そこで、上記線幅については固定値とせず、例えばユーザ設定パラメータとして設定変更が可能な状態としておくことが好ましい。
【0091】
図２０は探索範囲設定方法の第２の態様を示す図である。より具体的には、図２０（ａ）は探索範囲設定方法の第２の態様を示すフローチャートであり、図２０（ｂ）は処理過程における画像の例を模式的に示す図である。なお、概念の理解を容易にするために、図２０（ｂ）ではスフェアの外形として単純な楕円形が用いられている。この方法では、先に求められている輪郭情報に基づき、基準画像、比較画像のそれぞれにつき、スフェアの輪郭内部を表す画像Ｉ８，Ｉ９を描画する（ステップＳ６０１、Ｓ６０２）。そして、それらの画像Ｉ８，Ｉ９は排他的論理和合成される（ステップＳ６０３）。
【0092】
これにより得られる画像Ｉ１０では、２つの画像の一方でのみスフェアの輪郭内とされる領域のみが淡色（白）で表され、それ以外の領域は濃色（黒）で表される。この淡色領域を所定幅だけ膨張させる膨張処理を実行すると（ステップＳ６０４）、処理後の画像Ｉ１１には、基準画像および比較画像それぞれにおけるスフェアの輪郭（点線）を共に含み、かつそれらの輪郭の周囲まで広がる環状の帯状領域が現れる。この帯状領域が探索範囲Ｒｓとされる。ここでは画像Ｉ１１を第２の態様における「探索範囲画像」と称することとする。
【0093】
こうしていずれかの態様の探索範囲設定方法により生成される、探索範囲Ｒｓを表す探索範囲画像をマスクとして、収縮膨張画像Ｉ４ａに作用させる。このことは、収縮膨張画像Ｉ４ａと探索範囲画像Ｉ７または探索範囲画像Ｉ１１とを論理積合成することと等価である（ステップＳ４０８）。
【0094】
図２１は収縮膨張画像と探索範囲画像との論理積合成画像を示す図である。ここでは探索範囲画像として第１の態様の画像Ｉ７を用いて説明する。基準画像と比較画像との間におけるスフェアの外周部および内部の差分を表す収縮膨張画像Ｉ４ａと、探索範囲画像Ｉとを論理積合成することで、合成画像Ｉ１２では、スフェアの外周部の近傍における変化に対応する像のみが示され、収縮膨張画像Ｉ４ａに存在していたスフェア内部の像はマスクされて合成画像には現れない。これにより、外周部の位置変化、つまりスフェアの変位のみを選択的に抽出することが可能になる。この合成画像Ｉ１２を以下では「動き差分画像」と称する。
【0095】
変位の定量化は以下のようにして行われる。すなわち、動き差分画像Ｉ１２から、動き差分面積が求められる（ステップＳ４０９）。動き差分面積は、動き差分画像Ｉ１２における白色部分の面積として求められ、簡易的には該画像における白色画素の数により表すことができる。基準画像と比較画像との間におけるスフェアの変位が大きいほど動き差分面積は大きくなり、したがって動き差分面積は、変位の大きさを定量的に表す数値として利用可能である。
【0096】
ただし、動き差分面積はスフェアの大きさの影響を受ける。すなわち、面積が大きく外周部の周囲長の長いスフェアとより小さいスフェアとでは、たとえ変位量が同じであったとしても、変位の結果として現れる動き差分面積は、大きいスフェアの方が当然に大きくなる。これを是正し、スフェアの大きさに関わらず変位の大きさを表すことのできる指標として、以下の式：
Ｓ＝（動き差分面積）／（スフェアの平均輪郭長さ）
により、動きスコアＳを定義する。ここに、「スフェアの平均輪郭長さ」は、基準画像におけるスフェアの輪郭長さと比較画像におけるスフェアの輪郭長さとの平均値である。スフェアの輪郭長さについては、スフェア抽出処理の結果から導出することができ、例えば輪郭に相当する画素の数により表すことができる。
Ｓ＝（動き差分面積）／（スフェアの平均輪郭長さ）
により、動きスコアＳを定義する。ここに、「スフェアの平均輪郭長さ」は、基準画像におけるスフェアの輪郭長さと比較画像におけるスフェアの輪郭長さとの平均値である。スフェアの輪郭長さについては、スフェア抽出処理の結果から導出することができ、例えば輪郭に相当する画素の数により表すことができる。
【0097】
このように定義される動きスコアＳは、スフェアの変位により画像間に生じる差分がスフェアのサイズで正規化されたものであり、スフェアの大小に関わりなくその変位の大きさを指標するものとなる。
【0098】
基準画像と比較画像との一の組み合わせごとに、画像内の全てのスフェアにつき動きスコアＳが算出される（ステップＳ４１０）。そして、上記したように、基準画像と比較画像との組み合わせを順次変化させながら（ステップＳ４１１、Ｓ４１２）、それぞれの組み合わせでスフェアの動きスコアＳが算出される。このうち動きスコアＳが予め定められた閾値より大きいものについては、当該スフェアが有意な動きを示したものとみなすことができる。動きの検出されたスフェアについては、その動きの方向とともに記録されることが望ましい。
【0099】
図１に戻って、次に、動きスコアＳが算出された各スフェアにつき、その算出結果に基づき拍動判定が行われる（ステップＳ１０６）。スフェアの拍動の有無については、例えば以下のようにして判定することができる。
【0100】
図２２は拍動判定工程を示すフローチャートである。より具体的には、図２２は１つのスフェアについて拍動の有無を判定するための処理を示している。各スフェアには固有の識別番号を含むインデックスが付されているので、各フレームでの検出結果から同一識別番号のものを抽出することにより、１つのスフェアの動きを把握することができる。動き検出工程で求められた１つのスフェアについての動き検出結果を時系列順に配列すると（ステップＳ７０１）、当該スフェアが経時的にどのような動きを示したかが明らかになる。
【0101】
拍動はスフェアの継続的かつ周期的な運動である。したがって、スフェアが拍動していれば、当該スフェアについては複数回の動きが検出され、かつそれらの中には変位方向が反対のものが含まれているはずである。そこで、判定対象のスフェアがこのような条件を満たしていれば（ステップＳ７０２およびステップＳ７０３においてともにＹＥＳ）、当該スフェアは拍動していると判定される（ステップＳ７０４）。一方、動きの検出が１回のみである（ステップＳ７０２においてＮＯ）、あるいは複数回の動きの方向が同じである（ステップＳ７０３においてＮＯ）場合には、当該スフェアは拍動していないと判定される（ステップＳ７０５）。このようにして、各スフェアの拍動の有無を判定することができる。動き方向が反対であるか否かについては、例えば以下のようにして判断することができる。
【0102】
図２３は動き方向の判断例を示す図である。検出された動きの方向をベクトル表示したとき、図２３（ａ）に示すように、１つのベクトルＶ１の方向と、他の１つのベクトルＶ２の方向とのなす角θが９０度より小さければ、これらの動きは同方向のものと判断することができる。このような動きは拍動によるものではなく、単にスフェアが一方向へ移動しているものと考えられる。一方、図２３（ｂ）に示すように、２つのベクトルＶ１，Ｖ２の角θが９０度より大きいとき、これらの動きは反対方向のものと判断することができる。なおスフェアの動き方向については、上記した動き検出結果から求めてもよく、また例えばスフェアの重心位置の経時的な変化から求めてもよい。
【0103】
各スフェアの拍動に関する判定結果は、任意の態様で出力される（ステップＳ１０７）。例えば図１３（ｂ）に示すスフェア輪郭の画像Ｉｏにおいて、各スフェアの拍動の有無を識別できるような加工、例えば色分けを施して表示出力することが考えられる。また、図２に示すフレーム画像Ｉｆに対して、拍動するスフェア、拍動しないスフェアおよびスフェア以外のオブジェクトを互いに区別する加工を施してもよい。また、各スフェアの拍動の有無やその振幅、周期などの詳細な情報をインデックス情報とともにリストとして出力する態様であってもよい。その他、用途に応じて種々の態様で判定結果を出力することが可能である。
【0104】
以上のように、この実施形態の拍動検出処理では、撮像時刻が異なる基準画像と比較画像との画像内容の差分に基づいてスフェアの変位が求められ、その結果からスフェアの動きが検出される。拍動に起因する変位はスフェアの外周部の近傍に現れる一方、スフェア内部のテクスチャ変化は必ずしも拍動の有無を反映しない。このため、本実施形態では、抽出されたスフェアの輪郭の近傍に探索範囲を設定し、該範囲内での画像の差分を変位の大きさの指標としている。このため、内部テクスチャの変化に影響されることなく、拍動に起因するスフェアの変位を効率よく求めることが可能である。
【0105】
また、スフェアの拍動周期が必ずしも一定でなく、個体ばらつきもあることから、基準画像に組み合わせる比較画像については撮像時刻の異なる複数の画像を順次適用している。このため、種々の周期で生じ得るスフェアの拍動を確実に検出することが可能である。
【0106】
またスフェアの変位の大きさについては、変位によって差分画像に現れる領域の面積をスフェアの輪郭長さで正規化した動きスコアＳを用いて評価する。このため、スフェアの大きさの違いに対してロバストな検出および判定が可能となっている。
【0107】
＜装置構成例＞
図２４はこの実施形態の画像処理を実行可能なコンピュータ装置の構成例である。コンピュータ装置１は、例えばパーソナルコンピュータとして一般的な構成を有するものであり、ＣＰＵ（Central Processing Unit）１０、メモリ１４、ストレージ１５、入力デバイス１６、表示部１７、インターフェース１８およびディスクドライブ１９などを備えている。
図２４はこの実施形態の画像処理を実行可能なコンピュータ装置の構成例である。コンピュータ装置１は、例えばパーソナルコンピュータとして一般的な構成を有するものであり、ＣＰＵ（Central Processing Unit）１０、メモリ１４、ストレージ１５、入力デバイス１６、表示部１７、インターフェース１８およびディスクドライブ１９などを備えている。
【0108】
ＣＰＵ１０は、予め用意された制御プログラムを実行することで、上記した画像処理を実行するための機能ブロックとしての画像処理部１１をソフトウェア的に実現する。なお、画像処理部１１を実現するための専用ハードウェアが設けられてもよい。メモリ１４はＣＰＵ１０の演算過程で生成される各種データを一時的に記憶する。ストレージ１５は、ＣＰＵ１０が実行すべき制御プログラムのほか、原画像の画像データや処理後の画像データ等を長期的に記憶する。
【0109】
入力デバイス１６は、オペレータからの指示入力を受け付けるためのものであり、例えばマウス、キーボードなどを含む。また、表示部１７は画像を表示する機能を有する例えば液晶ディスプレイであり、原画像や処理後の画像、オペレータへのメッセージ等種々の情報を表示する。なお、入力デバイスと表示部とが一体化されたタッチパネルが設けられてもよい。
【0110】
インターフェース１８は、電気通信回線を介して外部装置との間で各種データ交換を行う。ディスクドライブ１９は、画像データや制御プログラム等各種のデータを記録した外部の記録ディスク２を受け入れる。記録ディスク２に記憶された画像データや制御プログラム等は、ディスクドライブ１９により読み出され、ストレージ１６に記憶される。ディスクドライブ１９はコンピュータ装置１内で生成されたデータを記録ディスク２に書き込む機能を備えていてもよい。
【0111】
本実施形態の画像処理をコンピュータ装置１に実行させるための制御プログラムについては、これを記録した記録ディスク２にディスクドライブ１９がアクセスして読み出される態様であってもよく、インターフェース１８を介して外部装置から与えられる態様であってもよい。原画像データについても同様である。
【0112】
＜その他＞
以上説明したように、本実施形態においては、スフェアＳｐが本発明の「心筋細胞のコロニー」および「コロニー領域」に相当している。また、連続撮像により取得される各フレーム画像Ｉｆを適宜の間隔で間引いた画像Ｉｆ１，Ｉｆ２，…が、本発明の「原画像」に相当している。拍動検出処理における基準画像Ｉ１および比較画像Ｉ２が、原画像の一部であるととともに、それぞれ本発明の「第１画像」および「第２画像」に相当する。
以上説明したように、本実施形態においては、スフェアＳｐが本発明の「心筋細胞のコロニー」および「コロニー領域」に相当している。また、連続撮像により取得される各フレーム画像Ｉｆを適宜の間隔で間引いた画像Ｉｆ１，Ｉｆ２，…が、本発明の「原画像」に相当している。拍動検出処理における基準画像Ｉ１および比較画像Ｉ２が、原画像の一部であるととともに、それぞれ本発明の「第１画像」および「第２画像」に相当する。
【0113】
なお、本発明は上記した実施形態に限定されるものではなく、その趣旨を逸脱しない限りにおいて上述したもの以外に種々の変更を行うことが可能である。例えば、上記実施形態では、試料の撮像、スフェアの抽出および拍動の判定までが一連の処理として記載されているが、常にこのような処理を要するものではない。例えば、過去に撮像された試料の画像を用いてスフェア抽出および拍動判定を行うようにしてもよい。また、スフェア抽出と拍動判定とについてもそれぞれ個別に実行されてもよい。
【0114】
また、上記実施形態では、画像内におけるスフェアの変位を求めることにより動き検出を行っているが、拍動によるスフェアの動きとしては周期的な収縮・弛緩として現れるものもある。上記実施形態は画像の差分の大きさから動きを検出するので、このような収縮・弛緩であっても検出可能である。
【0115】
また例えば、上記実施形態におけるスフェア抽出処理は、撮像された１つ１つの画像からスフェアのような比較的円形に近いオブジェクトを特定する処理であり、単独の静止画に対しても適用可能なものである。また、その対象物や利用目的も、心筋細胞スフェアやその拍動判定に限定されるものではなく、種々の細胞コロニーを対象として適用することが可能である。
【0116】
また例えば、上記実施形態における拍動判定処理を実行するのに際して、予め画像内のスフェアを特定しておく必要があるが、スフェアを特定するための処理については上記に限定されるものではなく、他の方法で特定されたスフェアの情報を用いて、上記の拍動判定処理を実行してもよい。
【0117】
またたとえば、上記実施形態では、撮像された各フレーム画像に対しスフェア抽出処理が実行され、フレームの間引きを行った上で拍動検出処理が実行される。これに代えて、例えば撮像された原画像に対しフレームの間引きを行い、間引き後のフレーム画像を用いてスフェア抽出処理および拍動検出処理が実行されてもよい。
【0118】
以上、具体的な実施形態を例示して説明してきたように、この発明に係る画像処理方法では、第１画像と第２画像との差分画像を２値化した２値化画像と、探索範囲に対応する探索範囲画像との論理積に相当する領域の面積に対応する値が、変位の大きさの指標値とされてもよい。このような構成によれば、コロニーの変位の大きさを差分として定量化することができ、しかも輪郭の近傍における変位のみを指標値に反映させることができるので、拍動を検出する目的に特に好適な指標値を得ることができる。
【0119】
この場合、第１画像における対象コロニー領域の輪郭長と第２画像における対象コロニー領域の輪郭長との平均値で面積を正規化した値が変位の大きさの指標値とされてもよい。このような構成によれば、変位量が同じでも面積の大きいコロニーほど差分が大きくなり指標値も大きくなるという問題を解消し、コロニーのサイズに対してロバストな判定が可能となる。
【0120】
また、複数の原画像から選択した一の第１画像に対し、撮像時刻が互いに異なる複数の原画像が１つずつ順に第２画像に設定され、それぞれの組み合わせについてコロニー領域の変位量が検出されてもよい。このような構成によれば、拍動周期が必ずしも一律でないスフェアについても、拍動の検出漏れを回避することができる。
【0121】
また、撮像時刻が互いに異なる複数の原画像が１つずつ順に第１画像に設定され、当該第１画像より後に撮像された原画像の一が第２画像として設定されてもよい。このような構成によれば、撮像中の任意のタイミングで生じるコロニー領域の動きを漏れなく検出することができる。
【0122】
また、第１画像と第２画像との間における対象コロニー領域の重心の移動方向を変位の方向とし、対象コロニー領域について、所定の閾値より大きい変位が複数回検出され、かつそれらの変位に、方向が互いに反対であるものが含まれる場合に、当該対象コロニー領域に対応するコロニーについて拍動ありと判定されるようにしてもよい。このような構成によれば、拍動とは異なる一方向への移動を誤って拍動と判定することが防止される。例えば２つの変位について、それぞれの変位方向を示す方向ベクトルのなす角が９０度より大きいときそれらの方向が反対であるとみなすことができる。
【0123】
また、原画像からコロニー領域を抽出する工程は、原画像から所定のエッジ強度を有する閉曲線を抽出する工程と、抽出された閉曲線で囲まれた領域を示す画像と原画像を２値化した画像との論理積に相当する合成画像を作成する工程と、合成画像に含まれるオブジェクトに対し分水嶺法によるオブジェクト分割を実行する工程とを含むように構成されてもよい。このような構成によれば、細胞のコロニーが高密度に含まれ互いに接触しているような試料においても、それらを的確に分離しつつ、それぞれのコロニーについて拍動の有無を個別に判定することが可能である。
【産業上の利用可能性】
【0124】
本発明は、細胞を撮像し画像を評価する生化学や医療の分野に適用可能であり、例えばｉＰＳ細胞やＥＳ細胞から分化培養される心筋細胞の評価に特に好適である。
【符号の説明】
【0125】
１ コンピュータ装置
２ 記録ディスク（記録媒体）
Ｉ１ 原画像（第１画像）
Ｉ２ 原画像（第２画像）
Ｉ３ 差分画像
Ｉ４ ２値化画像
Ｉ７，Ｉ１１ 探索範囲画像
Ｒｓ 探索範囲
Ｓｐ スフェア（コロニー、コロニー領域）
２ 記録ディスク（記録媒体）
Ｉ１ 原画像（第１画像）
Ｉ２ 原画像（第２画像）
Ｉ３ 差分画像
Ｉ４ ２値化画像
Ｉ７，Ｉ１１ 探索範囲画像
Ｒｓ 探索範囲
Ｓｐ スフェア（コロニー、コロニー領域）
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】