JP2020064349A - 画像処理方法、画像処理装置、および動作解析システム - Google Patents

Abstract

【課題】細胞集合体を同時に多数観察でき、且つ、細胞集合体の動く速さを精度良く観察できる技術を提供する。【解決手段】複数の細胞集合体を撮像したフレームが撮像された順番で並ぶ動画像を取得する工程と、前記動画像における前記細胞集合体の動作を解析する解析領域を、前記細胞集合体ごとに設定する工程と、前記解析領域を構成する画素ごとにｎ（ｎは１以上の自然数）番目のフレームにおける階調値とｎ＋ｍ（ｍは１以上の自然数）番目のフレームにおける階調値との差分を算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程と、前記細胞集合体の動く速さを表すパラメータとして前記差分の標準偏差および前記差分の絶対値平均から選ばれる少なくとも１つを前記細胞集合体ごとに算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程とを有する、画像処理方法。【選択図】図２

Description

本開示は、画像処理方法、画像処理装置、および動作解析システムに関する。

従来、細胞活動を測定する方法として、細胞活動によって生じる細胞周辺の電位変化を電極で測定する方法があった。この方法では、細胞を培養する培養容器として、電極が埋め込まれたものを用いる必要があった。

また、別の細胞活動を測定する方法として、細胞活動によって生じる細胞内部のイオン濃度変化を測定する方法があった。この方法では、例えばカルシウムイオンに結合することにより蛍光する蛍光色素を、細胞内部に加える必要があった。

特許文献１に記載の画像処理方法は、細胞が撮像された動画像を構成する複数のフレームのうち、現在時刻のフレームと、１つ前の時刻のフレームとの間での細胞の動き（位置の時間変化）を示す動きベクトルを算出する。この方法によれば、特殊な培養容器を使用することなく、また、細胞を傷付けることなく、細胞活動を観察できる。

特許文献２に記載の画像処理方法は、細胞が撮像された動画像を構成する複数のフレームのそれぞれと基準フレームとの差分画像を作成し、差分画像から例えば平均ピクセル強度を導出し、平均ピクセル強度から細胞の運動を定量化する。この方法によれば、特殊な培養容器を使用することなく、また、細胞を傷付けることなく、細胞活動を観察できる。

特開２０１２−１０５６３１号公報 特表２０１３−５３３９８６号公報

特許文献１によれば、細胞集合体を構成する複数個の細胞の動きベクトルを算出すべく、個々の細胞の輪郭を特定できる程度の高解像度の画像が必要であり、低倍率の画像で細胞集合体を同時に多数観察することが困難であった。

特許文献２によれば、基準フレームは、動画像を構成する複数のフレームから選択される。基準フレームが撮像された時刻と各フレームが撮像された時刻との時間差は、各フレームが撮像された時刻に応じて変動する。その時間差が長過ぎると、細胞の動く方向が途中で反転する。それゆえ、差分画像から細胞の動きがほとんど無いと解析された場合に、実際には細胞の動きが活発なことがあった。つまり、特許文献２によれば、差分画像から、細胞の動く速さを精度良く観察できなかった。

本開示の一態様は、細胞集合体を同時に多数観察でき、且つ、細胞集合体の動く速さを精度良く観察できる、技術を提供する。

本開示の一態様に係る画像処理方法は、
複数の細胞集合体を撮像したフレームが撮像された順番で並ぶ動画像を取得する工程と、
前記動画像における前記細胞集合体の動作を解析する解析領域を、前記細胞集合体ごとに設定する工程と、
前記解析領域を構成する画素ごとにｎ（ｎは１以上の自然数）番目のフレームにおける階調値とｎ＋ｍ（ｍは１以上の自然数）番目のフレームにおける階調値との差分を算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程と、
前記細胞集合体の動く速さを表すパラメータとして前記差分の標準偏差および前記差分の絶対値平均から選ばれる少なくとも１つを前記細胞集合体ごとに算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程とを有する。

本開示の一態様によれば、細胞集合体を同時に多数観察でき、且つ、細胞集合体の動く速さを精度良く観察できる。

図１は、一実施形態に係る動作解析システムを示す図である。 図２は、収縮完了時および弛緩開始時の細胞集合体を通過する光と、その光を受光する受光素子との関係の一例を示す図である。 図３は、弛緩完了時および収縮開始時の細胞集合体を通過する光と、その光を受光する受光素子との関係の一例を示す図である。 図４は、撮像装置によって撮像される動画像と、動画像を構成する複数のフレームとの関係を示す図である。 図５は、一実施形態に係る画像処理装置の構成要素を機能ブロックで示す図である。 図６は、ｎ番目のフレームの一例を示す図である。 図７は、受光素子と、受光素子で受光される光の強度に相当する階調値との関係の一例を示す図である。 図８は、一実施形態に係る画像処理方法を示すフローチャートである。 図９は、試験例１に係るＢ（ｎ）の時間変化を示すグラフである。 図１０は、試験例１に係るＳ（ｎ）の時間変化を示すグラフである。 図１１は、試験例１に係るΔＢ（ｎ）の時間変化を示すグラフである。 図１２は、試験例１に係るΔＳ（ｎ）の時間変化を示すグラフである。 図１３は、試験例２に係るΔＳ（ｎ）とΔＡ（ｎ）の時間変化を示すグラフである。 図１４は、試験例３に係るΔＳ（ｎ）の時間変化を示すグラフである。 図１５は、試験例３に係るΔＰｘｙの平均値の時間変化を示すグラフであって、基準フレームとして９５番目のフレームを選択した時のグラフである。 図１６は、試験例３に係るΔＰｘｙの平均値の時間変化を示すグラフであって、基準フレームとして１番目のフレームを選択した時のグラフである。

以下、本開示の実施形態について図面を参照して説明する。なお、各図面において同一の又は対応する構成には同一の又は対応する符号を付し、説明を省略することがある。

図１は、一実施形態に係る動作解析システムを示す図である。動作解析システム１０は、複数の細胞が集合した細胞集合体の動作を解析する。動作解析システム１０は、複数の細胞集合体の動画像を撮像する撮像装置２０と、細胞集合体の動作を解析すべく、撮像装置２０によって撮像された動画像を処理する画像処理装置３０とを備える。動作解析システム１０は、画像処理装置３０の解析結果を表示する画像表示装置４０をさらに備える。

撮像装置２０は、複数の細胞集合体の動画像を撮像する。複数の細胞集合体は、細胞集合体の損傷を抑制する目的で、培養容器１２に収容された状態で撮像される。培養容器１２は、細胞集合体を培養するものである。培養容器１２は、一般的なものであるので、説明を省略する。撮像素子２１が細胞集合体の動画像を撮像する間、培養容器１２と撮像素子２１との位置関係は固定される。

細胞集合体は、例えば、スフェロイド（Spheroid）と呼ばれる、三次元的な細胞のコロニーである。三次元的な細胞のコロニーは、二次元的な細胞のコロニーに比べて、生体に類似した構造を有する。それゆえ、創薬スクリーニングなどへの利用が期待されている。

なお、細胞集合体は、本実施形態では三次元的な細胞のコロニーであるが、二次元的な細胞のコロニーであってもよい。

撮像装置２０は、複数の細胞集合体を撮像する撮像素子２１を有する。撮像素子２１は、例えばＣＭＯＳイメージセンサまたはＣＣＤイメージセンサである。撮像素子２１は、細胞集合体を通過した光を受光することにより、受光した光の強度に応じた量の電荷を蓄積する受光素子２２を複数有する。

複数の受光素子２２は、行列状に配置される。受光素子２２は、撮像素子２１の画素ごとに１つずつ配置される。受光素子２２で受光される光の強度が強いほど、受光素子２２に蓄積される電荷量の大きさが大きい。

図２は、収縮完了時および弛緩開始時の細胞集合体を通過する光と、その光を受光する受光素子との関係の一例を示す図である。図３は、弛緩完了時および収縮開始時の細胞集合体を通過する光と、その光を受光する受光素子との関係の一例を示す図である。図２および図３において、矢印は光の進行方向を表す。

細胞集合体１４は、例えば数百個の細胞１５で構成される。細胞集合体１４は、例えば主に心筋細胞を含み、図２および図３に示すように収縮と弛緩とを繰り返す。収縮と弛緩とを繰り返すことを、拍動とも呼ぶ。

細胞集合体１４に含まれる細胞１５が動く時、細胞１５の重なり具合が変化し、細胞集合体１４を通過する光の伝播方向が変化する。細胞１５の重なり具合によって、光が強く回折される場所と、光がほとんど回折されずに直進する場所とが生じるからである。

なお、細胞集合体１４は、心筋細胞を含むものには限定されない。細胞集合体１４は、動きのある細胞１５を含むものであればよい。細胞集合体１４に含まれる細胞１５が動けば、細胞１５の重なり具合が変化するからである。

細胞集合体１４に含まれる細胞１５が動く時、細胞１５の重なり具合が変化する。その結果、細胞集合体１４を通過する光の方向が変化するので、受光素子２２で受光される光の強度が変化し、受光素子２２に蓄積される電荷量が変化する。

図４は、撮像装置によって撮像される動画像と、動画像を構成する複数のフレームとの関係を示す図である。図４に示すように、動画像５０は、撮像された順番で並ぶ複数のフレーム５１を有する。フレーム５１とは、静止画像のことである。

フレーム５１は、行列状に並ぶ複数の画素５２を有する。フレーム５１の画素５２と、撮像素子２１の画素とは一対一で対応する。撮像素子２１の画素は、受光素子２２である。

複数の画素５２の階調値は、互いに異なる受光素子２２に蓄積される電荷量で決められる。つまり、複数の画素５２の階調値は、互いに異なる受光素子２２で受光される光の強度で決められる。

フレーム５１は、例えば白黒の濃淡で表現されるモノクロ画像であり、複数の受光素子２２で受光される光の強度を表す。受光素子２２で受光される光の強度が強くなるほど、画素５２の色が黒色から白色に変化する。

画素５２の階調は、特に限定されないが、例えば８ビット（２５６階調）である。この場合、画素５２の階調値は、０〜２５５の整数である。受光素子２２で受光される光の強度が強いほど、画素５２の階調値が大きくなる。階調値が０である画素５２の色は黒色であり、階調値が２５５である画素５２の色は白色である。

フレームレートは、ｎ（ｎは１以上の自然数）番目のフレーム５１（ｎ）とｎ＋１番目のフレーム５１（ｎ＋１）との間隔が細胞集合体１４の動作の周期（例えば拍動の周期）よりも十分に短くなるように設定される。フレームレートは、単位時間（例えば１秒）当たりのフレーム５１の数のことである。

動画像５０に写る細胞集合体１４の個数は、撮像装置２０の撮影倍率などで決まる。撮影倍率とは、撮像素子２１上での被写体の長さを、実際の被写体の長さで除した値である。撮影倍率が小さいほど、同時に観察できる細胞集合体１４の個数が増える反面、細胞集合体１４を構成する個々の細胞１５が不鮮明になる。

撮像装置２０の撮影倍率は、例えば４倍未満であり、好ましくは２倍以下である。撮像装置２０の撮影倍率は、好ましくは０．３倍以上、より好ましくは０．５倍以上である。なお、同時に観察できる細胞集合体１４の個数は、撮像装置２０の撮影倍率の他、撮像素子２１の面積で決まる。

細胞集合体１４に照射される光の波長は、細胞集合体１４の損傷を抑制するように設定され、例えば６２０ｎｍ〜７５０ｎｍに設定される。また、細胞集合体１４に照射される光の強度は、細胞集合体１４の損傷を抑制するように設定され、例えば１．０×１０^−１４Ｗ／μｍ^２〜１．２×１０^−８Ｗ／μｍ^２に設定される。

画像処理装置３０は、例えばコンピュータで構成され、図１に示すように、ＣＰＵ（Central Processing Unit）３１と、メモリなどの記憶媒体３２とを備える。記憶媒体３２には、動作解析システム１０において実行される各種の処理を制御するプログラムが格納される。画像処理装置３０は、記憶媒体３２に記憶されたプログラムをＣＰＵ３１に実行させることにより、動作解析システム１０の動作を制御する。また、画像処理装置３０は、入力インターフェース３３と、出力インターフェース３４とを備える。画像処理装置３０は、入力インターフェース３３で外部からの信号を受信し、出力インターフェース３４で外部に信号を送信する。

図５は、一実施形態に係る画像処理装置の構成要素を機能ブロックで示す図である。図５に図示される各機能ブロックは概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。各機能ブロックの全部または一部を、任意の単位で機能的または物理的に分散・統合して構成することが可能である。各機能ブロックにて行われる各処理機能は、その全部または任意の一部が、ＣＰＵ３１にて実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現されうる。

画像処理装置３０は、撮像装置２０によって撮像された動画像５０を取得する画像取得部３５を有する。画像取得部３５は、有線または無線で、撮像装置２０によって撮像された動画像５０を取得する。動画像５０は撮像された順番で並ぶ複数のフレーム５１を有し、複数のフレーム５１のそれぞれは複数の細胞集合体１４を撮像したものである。複数のフレーム５１は、それぞれ、撮像された順番と対応付けて記憶媒体３２に記憶される。

図６は、ｎ番目のフレームの一例を示す図である。図６には、ｎ番目のフレーム５１（ｎ）の一部のみを示す。ｎ番目のフレーム５１（ｎ）に写る細胞集合体１４の総数は、図６に示す細胞集合体１４の数よりも多い。図６において、白色の破線は、細胞集合体１４の動作を解析する解析領域５３を表す。

画像処理装置３０は、動画像５０における細胞集合体１４の動作を解析する解析領域５３を、細胞集合体１４ごとに設定する解析領域設定部３６を有する。解析領域設定部３６は、細胞集合体１４の輪郭の内部に解析領域５３を設定する。解析領域５３は、時間の経過に関係なく、常に、細胞集合体１４の輪郭の内部に収まるように設定される。

解析領域設定部３６は、例えばユーザによる入力操作を受け付ける入力装置からの信号に基づいて解析領域５３を設定する。入力装置は、ユーザによる入力操作に応じた信号を、解析領域設定部３６に送信する。なお、解析領域設定部３６は、パターンマッチングまたは二値化処理などの画像処理によって細胞集合体１４の輪郭を特定し、細胞集合体１４の輪郭の内部に解析領域５３を設定してもよい。

解析領域設定部３６は、例えば、細胞集合体１４ごとに、解析領域５３を１つずつ設定する。解析領域５３の大きさは特に限定されないが、例えば３０×３０個の画素５２で解析領域５３が構成される。解析領域５３の形状は、本実施形態では正方形であるが、例えば長方形、三角形、円形、または楕円形などでもよく、特に限定されない。

なお、本実施形態では１つの細胞集合体１４の輪郭の内部に１つの解析領域５３が設定されるが、本開示の技術はこれに限定されない。１つの細胞集合体１４の輪郭の内部に、複数の解析領域５３が間隔をおいて設定されてもよい。

画像処理装置３０は、解析領域５３を構成する画素５２ごとに、ｎ番目のフレーム５１（ｎ）における階調値Ｇ（ｎ）とｎ＋ｍ番目のフレーム５１（ｎ＋ｍ）における階調値Ｇ（ｎ＋ｍ）との差分ΔＧ（ｎ）を算出する差分算出部３７を有する。ｎは１以上の自然数であり、ｍは１以上の自然数である。

ｎ番目のフレーム５１（ｎ）が撮像された時刻Ｔ（ｎ）は、ｎ／Ｆで表される。また、ｎ＋ｍ番目のフレーム５１（ｎ＋ｍ）が撮像された時刻Ｔ（ｎ＋ｍ）は、（ｎ＋ｍ）／Ｆで表される。ここで、Ｆは、フレームレートである。

１つの解析領域５３がｉ×ｊ個の画素５２で構成される場合、解析領域５３ごとにｉ×ｊ個のΔＧ（ｎ）が算出される。ｉは、１以上の自然数であり、例えば３０である。ｊは、１以上の自然数であり、例えば３０である。ｉ×ｊは、２以上の自然数であり、例えば９００である。

ΔＧ（ｎ）は、解析領域５３の内部の画素５２について算出される。ΔＧ（ｎ）は、解析領域５３の外部の画素５２について、算出されてもよいが、算出されなくてよい。無駄な演算が行われないので、ＣＰＵ３１の演算量を低減できる。

ΔＧ（ｎ）は、「ΔＧ（ｎ）＝Ｇ（ｎ）−Ｇ（ｎ＋ｍ）」の式から求められる値である。なお、ΔＧ（ｎ）は、「ΔＧ（ｎ）＝Ｇ（ｎ＋ｍ）−Ｇ（ｎ）」の式から求められる値でもよい。ΔＧ（ｎ）が正であっても負であっても、ΔＧ（ｎ）の絶対値が変わらなければ、ΔＧ（ｎ）の標準偏差ΔＳ（ｎ）およびΔＧ（ｎ）の絶対値平均ΔＡ（ｎ）も変わらないからである。

差分算出部３７は、ΔＧ（ｎ）の算出を、ｍを固定すると共にｎを変更しながら、繰り返し行う。ΔＧ（ｎ）の時間変化を求めることができる。ΔＧ（ｎ）の時刻Ｔ（ｎ）は、ｎ／Ｆで表される。

例えば、差分算出部３７は、ΔＧ（ｎ）の算出を、ｍを固定すると共にｎを変更しながら繰り返し行うことにより、ｎが１つずつ異なるΔＧ（ｎ）の集合を得る。ｎが例えば２つずつ異なるΔＧ（ｎ）の集合を得る場合に比べて、ΔＧ（ｎ）の時間変化を詳細に求めることができる。なお、ｍは、詳しくは後述するが、好ましくは２以上である。

画像処理装置３０は、細胞集合体１４の動く速さを表すパラメータとして、ΔＧ（ｎ）の標準偏差ΔＳ（ｎ）およびΔＧ（ｎ）の絶対値平均ΔＡ（ｎ）から選ばれる少なくとも１つを、細胞集合体１４ごとに算出する速度算出部３８を有する。細胞集合体１４の動く速さを表すパラメータとして、ΔＳ（ｎ）およびΔＡ（ｎ）に着目した理由について図７を参照して説明する。

図７は、受光素子と、受光素子で受光される光の強度に相当する階調値との関係の一例を示す図である。図７（ａ）は、受光素子とＧ（ｎ）との関係の一例を示す図である。図７（ｂ）は、受光素子とＧ（ｎ＋ｍ）との関係の一例を示す図である。図７（ｃ）は、受光素子とΔＧとの関係の一例を示す図である。

図７に示す撮像素子２１は、７つの受光素子２２−１〜２２−７を有する。７つの受光素子２２−１〜２２−７は、左右方向に等ピッチで一列に並ぶ。そのピッチはＬである。７つの受光素子２２−１〜２２−７は、細胞集合体１４を通過する光を受光する。

図７（ａ）に示すように、時刻Ｔ（ｎ）において、１つの受光素子２２−２のみが細胞集合体１４を通過する光を受光し、残りの６つの受光素子２２−１、２２−３〜２２−７は細胞集合体１４を通過する光を受光しない。時刻Ｔ（ｎ）において受光素子２２−２で受光される光の強度は、画素５２の階調値Ｇ_０に相当する強度である。

時刻Ｔ（ｎ）において１つの受光素子２２−２で受光する光は、時刻Ｔ（ｎ）から時刻Ｔ（ｎ＋ｍ）までの時間ΔＴ（ΔＴ＝ｍ／Ｆ）の間に右方向に速さｖ（例えば、ｖ＜Ｌ／ΔＴ＝Ｌ×Ｆ／ｍ）で移動する。その結果、時刻Ｔ（ｎ）において１つの受光素子２２−２で受光する光は、時刻Ｔ（ｎ＋ｍ）において元の受光素子２２−２とその右隣の受光素子２２−３との両方にまたがって受光される。受光位置が移動するのは、図２および図３に示すように細胞集合体１４が動くからである。細胞集合体１４の動く速さが大きいほど、受光位置の移動する速さｖが大きい。

図７（ｂ）に示すように、時刻Ｔ（ｎ＋ｍ）において、２つの受光素子２２−２、２２−３のみが細胞集合体１４を通過する光を受光し、残りの５つの受光素子２２−１、２２−４〜２２−７は細胞集合体１４を通過する光を受光しない。時刻Ｔ（ｎ＋ｍ）において、１つの受光素子２２−２で受光される光の強度は、画素５２の階調値Ｇ_１に相当する強度である。また、時刻Ｔ（ｎ＋ｍ）において、別の１つの受光素子２２−３で受光される光の強度は、画素５２の階調値Ｇ_２に相当する強度である。

時刻Ｔ（ｎ＋ｍ）において２つの受光素子２２−２、２２−３で受光される光の強度の和は、時刻Ｔ（ｎ）において１つの受光素子２２−２で受光される光の強度に等しい。従って、Ｇ_１とＧ_２の和は、Ｇ_０に等しい。ΔＴが一定の場合、ｖが大きいほど、Ｇ_１が小さく、Ｇ_２が大きい。Ｇ_１とＧ_２との関係は、図７（ｂ）に示す式で表される。

図７（ａ）および図７（ｂ）に示すように、時刻Ｔ（ｎ）から時刻Ｔ（ｎ＋ｍ）までのΔＴの間に、２つの受光素子２２−２、２２−３のみにおいて受光量が変化し、残りの５つの受光素子２２−１、２２−４〜２２−７において受光量が変化しない。従って、図７（ｃ）に示すように、２つのΔＧ（ｎ）のみが０以外の数値ΔＧ_１、ΔＧ_２になり、残りの５つのΔＧ（ｎ）は０になる。ΔＧ_１とΔＧ_２の和は、０になる。ΔＴが一定の場合、ｖが大きいほど、ΔＧ_１が大きく、ΔＧ_２が小さい。ΔＧ_１とΔＧ_２との関係は、図７（ｃ）に示す式で表される。

図７（ｃ）に示すΔＧ（ｎ）の標準偏差ΔＳ（ｎ）は、下記式（１）で表される。
ΔＳ（ｎ）＝Ｇ_０×（ｖ×ΔＴ）／Ｌ×（２／７）^０．５・・・（１）
上記式（１）から明らかなように、ΔＳ（ｎ）は、受光位置の移動の速さｖで決まるので、細胞集合体１４の動きの速さで決まる。細胞集合体１４の動く速さが大きいほど、受光位置の移動の速さｖが大きく、ΔＳ（ｎ）が大きい。従って、ΔＳ（ｎ）を、細胞集合体１４の動く速さを表すパラメータとして使用できる。

図７（ｃ）に示すΔＧ（ｎ）の絶対値平均ΔＡ（ｎ）は、下記式（２）で表される。
ΔＡ（ｎ）＝Ｇ_０×（ｖ×ΔＴ）／Ｌ×（２／７）・・・（２）
上記式（２）から明らかなように、ΔＡ（ｎ）は、受光位置の移動の速さｖで決まるので、細胞集合体１４の動きの速さで決まる。細胞集合体１４の動く速さが大きいほど、受光位置の移動の速さｖが大きく、ΔＡ（ｎ）が大きい。従って、ΔＡ（ｎ）を、細胞集合体１４の動く速さを表すパラメータとして使用できる。

一方、図７（ｃ）に示すΔＧ（ｎ）の相加平均ΔＢ（ｎ）は、下記式（３）で表される。
ΔＢ（ｎ）＝０・・・（３）
上記式（３）から明らかなように、ΔＢ（ｎ）は、受光位置の移動の速さｖには依存しないので、細胞集合体１４の動きの速さには依存しない。従って、ΔＢ（ｎ）は、細胞集合体１４の動く速さを表すパラメータとして不適である。

以上説明したように、ΔＧ（ｎ）の統計値として、一部の統計値のみが、細胞集合体１４の動く速さを表すパラメータとして使用できる。

本実施形態の速度算出部３８は、ΔＧ（ｎ）の標準偏差ΔＳ（ｎ）およびΔＧ（ｎ）の絶対値平均ΔＡ（ｎ）から選ばれる少なくとも１つを算出する。ΔＳ（ｎ）およびΔＡ（ｎ）は細胞集合体１４の動く速さで決まるので、細胞集合体１４の動く速さを精度良く観察することができる。

ところで、上記特許文献２では、差分画像の平均ピクセル強度を算出する。上記特許文献２では、差分画像の作成に用いる基準フレームが固定されるので、差分画像は速さに関するデータを失ってしまうことがある。基準フレームが撮像された時刻と各フレームが撮像された時刻との時間差は、各フレームが撮像された時刻に応じて変動する。その時間差が長過ぎると、細胞１５の動く方向が途中で反転するからである。上記特許文献２の差分画像は速さに関するデータを失うことがあるので、その差分画像の平均ピクセル強度も当然に速さに関するデータを失うことがある。

これに対し、本実施形態では、ΔＧ（ｎ）の標準偏差ΔＳ（ｎ）およびΔＧ（ｎ）の絶対値平均ΔＡ（ｎ）から選ばれる少なくとも１つを算出する。ΔＧ（ｎ）は、Ｇ（ｎ）とＧ（ｎ＋ｍ）との差分であり、解析領域５３を構成する画素５２ごとに算出される。ΔＧ（ｎ）の算出は、ｍを固定すると共にｎを変更しながら、繰り返し行われる。本実施形態では、ｎ番目のフレーム５１が撮像された時刻Ｔ（ｎ）から、ｎ＋ｍ番目のフレーム５１が撮像された時刻Ｔ（ｎ＋ｍ）までの時間ΔＴは、ｎに関係なく、ｍ／Ｆである。ｍおよびＦは定数であるので、ΔＴは定数である。ΔＴがｎに応じて変動しないので、ｎの変動が原因でΔＧ（ｎ）が速さに関するデータを失うことはない。従って、本実施形態によれば、細胞集合体１４の動く速さを精度良く観察できる。また、本実施形態によれば、上記特許文献１の下記の問題をも解決できる。

上記特許文献１では、細胞集合体１４を構成する複数個の細胞１５の動きベクトルを算出すべく、個々の細胞１５の輪郭を特定できる程度の高解像度の画像が必要であり、低倍率の画像で細胞集合体１４を同時に多数観察することが困難であった。

これに対し、本実施形態では、細胞集合体１４の動きの速さを表すパラメータとしてΔＳ（ｎ）およびΔＡ（ｎ）から選ばれる少なくとも１つを算出する。ΔＳ（ｎ）およびΔＡ（ｎ）はΔＧ（ｎ）のバラツキを表す量であるので、ΔＧ（ｎ）の数は好ましくは１００以上である。ΔＧ（ｎ）の数が細胞集合体１４を構成する細胞１５の数と同程度で足りるので、個々の細胞１５の輪郭を特定できる程度の高解像度の画像は不要である。従って、本実施形態によれば、低倍率の画像で細胞集合体１４を同時に多数観察することができる。

１つの細胞集合体１４の動きの速さを表すΔＳ（ｎ）およびΔＡ（ｎ）は、その細胞集合体１４に設定される１つ以上の解析領域５３に含まれる画素５２の数と同数のΔＧ（ｎ）から求められる。１つの細胞集合体１４の動きの速さを表すΔＳ（ｎ）およびΔＡ（ｎ）の算出に用いられるΔＧ（ｎ）の数は、例えば１００以上であり、好ましくは１０００以上である。

本実施形態の速度算出部３８は、ΔＳ（ｎ）およびΔＡ（ｎ）から選ばれる少なくとも１つの算出を、ｍを固定すると共にｎを変更しながら、繰り返し行う。細胞集合体１４の動く速さの時間変化を精度良く観察することができる。また、本実施形態によれば、上記特許文献２の下記の問題をも解決できる。

上記特許文献２では、差分画像の作成に用いる基準フレームを、時系列に並ぶ複数のフレームから選択する必要がある。選択した基準フレームの撮像された時刻が図２に示す収縮完了時の時刻か図３に示す弛緩完了時の時刻かによって、差分画像が全く変わってしまう。従って、基準フレームの選択の仕方で、差分画像の平均ピクセル強度の波形が全く変わってしまう。

これに対し、本実施形態では、基準フレームを選択する必要がない。ｎ番目のフレーム５１が撮像された時刻Ｔ（ｎ）から、ｎ＋ｍ番目のフレーム５１が撮像された時刻Ｔ（ｎ＋ｍ）までの時間ΔＴは、ｎに関係なく、ｍ／Ｆである。ｍおよびＦは定数であるので、ΔＴは定数である。ΔＴが変動しないので、ｎの変動が原因でΔＳ（ｎ）の波形およびΔＡ（ｎ）の波形が変わることはない。

ところで、細胞集合体１４の動きを滑らかに撮像するには、単位時間当たりのフレーム５１の数（つまり、フレームレートＦ）を大きくすればよい。フレームレートＦが大きいと、時刻Ｔ（ｎ）から時刻Ｔ（ｎ＋１）までの時間（１／Ｆ）が短いので、その時間（１／Ｆ）の間の細胞集合体１４の動きの大きさが小さく、その動きの検出精度が低下する。

そこで、本実施形態の速度算出部３８は、差分ΔＧ（ｎ）を算出することを、ｍを２以上の自然数に固定すると共にｎを変更しながら繰り返し行うことにより、ｎが１つずつ異なるΔＧ（ｎ）の集合を得る。ｍを２以上の自然数に固定することにより、フレームレートＦの増加に伴うΔＴ（ΔＴ＝ｍ／Ｆ）の短縮を抑制できるので、ΔＴの間の細胞集合体１４の動きの大きさの低下を抑制できる。従って、細胞集合体１４の動きを滑らかに撮像することと、ΔＴの間の細胞集合体１４の動きの速さを精度良く検出することとを、両立できる。

ｍは、好ましくはｍ／Ｆの値が０．０１秒〜０．１秒を満たすように設定される。また、ｍは、例えば「ｍ＜Ｌ×Ｆ／ｖ」の式を満たすように設定される。「ｍ＜Ｌ×Ｆ／ｖ」の式が成立することは、時刻Ｔ（ｎ）において１つの受光素子２２−２で受光する光が時刻Ｔ（ｎ＋ｍ）において元の受光素子２２−２とその隣の受光素子２２−３との両方にまたがって受光されることを意味する（図７参照）。

画像処理装置３０は、速度算出部３８の算出結果を、画像表示装置４０に表示する表示処理部３９を有する。ユーザは、画像表示装置４０に表示された速度算出部３８の算出結果を見ることにより、細胞集合体１４の動きの速さを知ることができる。画像表示装置４０としては、フラットパネルディスプレイなどが用いられる。

図８は、一実施形態に係る画像処理方法を示すフローチャートである。図８に示す各工程は、画像処理装置３０によって実施される。なお、図８に示す工程の順序は、特に限定されない。例えば、工程Ｓ１０１〜Ｓ１０４は、並行して行われてもよい。

画像処理方法は、複数の細胞集合体１４を撮像したフレーム５１が撮像された順番で並ぶ動画像５０を取得する工程Ｓ１０１を有する。この工程Ｓ１０１では、画像取得部３５が、撮像装置２０によって撮像された動画像５０を取得する。動画像５０は複数のフレーム５１で構成され、複数のフレーム５１はそれぞれ撮像された順番と対応付けて記憶媒体３２に記憶される。

画像処理方法は、動画像５０における細胞集合体１４の動作を解析する解析領域５３を、細胞集合体１４ごとに設定する工程Ｓ１０２を有する。この工程Ｓ１０２では、解析領域設定部３６が、解析領域５３を細胞集合体１４ごとに設定する。解析領域設定部３６は、細胞集合体１４の輪郭の内部に解析領域５３を設定する。

画像処理方法は、解析領域５３を構成する画素ごとに階調値の差分ΔＧ（ｎ）を算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程Ｓ１０３を有する。この工程Ｓ１０３では、差分算出部３７が、ΔＧ（ｎ）の算出を、ｍを固定すると共にｎを変更しながら、繰り返し行う。例えば、差分算出部３７は、ΔＧ（ｎ）の算出を、ｍを固定すると共にｎを変更しながら繰り返し行うことにより、ｎが１つずつ異なるΔＧ（ｎ）の集合を得る。なお、ｍは、好ましくは２以上である。

画像処理方法は、細胞集合体１４の動く速さを表すパラメータとしてΔＳ（ｎ）およびΔＡ（ｎ）から選ばれる少なくとも１つを細胞集合体１４ごとに算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程Ｓ１０４を有する。この工程Ｓ１０４では、速度算出部３８が、ΔＳ（ｎ）およびΔＡ（ｎ）から選ばれる少なくとも１つを細胞集合体１４ごとに算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う。

ΔＳ（ｎ）およびΔＡ（ｎ）は、上記式（１）および上記式（２）から明らかなように細胞集合体１４の動く速さで決まる。また、ΔＳ（ｎ）およびΔＡ（ｎ）はΔＧ（ｎ）のバラツキを表す量であり、ΔＧ（ｎ）はＧ（ｎ）とＧ（ｎ＋ｍ）との差分である。時刻Ｔ（ｎ）から時刻Ｔ（ｎ＋ｍ）までの時間ΔＴは、ｎに関係なく、ｍ／Ｆである。ｍおよびＦは定数であるので、ΔＴは定数である。ΔＴがｎに応じて変動しないので、ｎの変動が原因でΔＧ（ｎ）が速さに関するデータを失うことはない。従って、本実施形態によれば、細胞集合体１４の動く速さを精度良く観察できる。

また、ΔＳ（ｎ）およびΔＡ（ｎ）はΔＧ（ｎ）のバラツキを表す量であるので、ΔＧ（ｎ）の数は好ましくは１００以上である。ΔＧ（ｎ）の数が細胞集合体１４を構成する細胞１５の数と同程度で足りるので、個々の細胞１５の輪郭を特定できる程度の高解像度の画像は不要である。従って、本実施形態によれば、低倍率の画像で細胞集合体１４を同時に多数観察することができる。

速度算出部３８は、例えば、差分ΔＧ（ｎ）を算出することを、ｍを２以上の自然数に固定すると共にｎを変更しながら繰り返し行うことにより、ｎが１つずつ異なるΔＧ（ｎ）の集合を得る。ｍを２以上の自然数に固定することにより、フレームレートＦの増加に伴うΔＴ（ΔＴ＝ｍ／Ｆ）の短縮を抑制できるので、ΔＴの間の細胞集合体１４の動きの大きさの低下を抑制できる。従って、細胞集合体１４の動きを滑らかに撮像することと、ΔＴの間の細胞集合体１４の動きの速さを精度良く検出することとを、両立できる。

画像処理方法は、上記工程Ｓ１０４の算出結果を、画像表示装置４０に表示する工程Ｓ１０５を有する。ユーザは、画像表示装置４０に表示された上記工程Ｓ１０４の算出結果を見ることにより、細胞集合体１４の動きの速さを知ることができる。

以下、実際に計測した統計値について、図９〜図１６を参照して説明する。

試験例１では、約２００個の細胞１５で構成される細胞集合体１４を１００個同時に撮像した動画像５０からＧ（ｎ）およびΔＧ（ｎ）を算出し、細胞集合体１４ごとにＢ（ｎ）、Ｓ（ｎ）、ΔＢ（ｎ）およびΔＳ（ｎ）を算出した。ここで、Ｂ（ｎ）はＧ（ｎ）の相加平均であり、Ｓ（ｎ）はＧ（ｎ）の標準偏差である。また、ΔＢ（ｎ）はΔＧ（ｎ）の相加平均であり、ΔＳ（ｎ）はΔＧ（ｎ）の標準偏差である。９６個の細胞集合体１４の代表として、３個の細胞集合体１４の結果を、図９〜図１２に実線と破線と一点鎖線とで示す。

図９は、試験例１に係るＢ（ｎ）の時間変化を示すグラフである。図１０は、試験例１に係るＳ（ｎ）の時間変化を示すグラフである。図１１は、試験例１に係るΔＢ（ｎ）の時間変化を示すグラフである。図１２は、試験例１に係るΔＳ（ｎ）の時間変化を示すグラフである。

図９から明らかなように、３個のＢ（ｎ）の波形は特徴的なパターンを有さないので、Ｂ（ｎ）から細胞集合体１４の動きを把握できないことが分かる。また、図１０から明らかなように、３個のＳ（ｎ）の波形は特徴的なパターンを有さないので、Ｓ（ｎ）から細胞集合体１４の動きを把握できないことが分かる。さらに、図１１から明らかなように、３個のΔＢ（ｎ）の波形は特徴的なパターンを有さないので、ΔＢ（ｎ）から細胞集合体１４の動きを把握できないことが分かる。

一方、図１２から明らかなように、３個のΔＳ（ｎ）の波形は、いずれも、特徴的なパターンを有し、具体的には大きいピークと小さいピークとの組合せを繰り返し有する。それゆえ、ΔＳ（ｎ）から細胞集合体１４の動きを把握できることが分かる。図１２において隣り合う２つのピークのうち、大きいピークは細胞集合体１４が収縮中であることを表し、小さいピークは細胞集合体１４が弛緩中であることを表す。

試験例２では、約２００個の細胞１５で構成される細胞集合体１４を１００個同時に撮像した動画像５０からＧ（ｎ）およびΔＧ（ｎ）を算出し、細胞集合体１４ごとにΔＳ（ｎ）およびΔＡ（ｎ）を算出した。ΔＳ（ｎ）はΔＧ（ｎ）の標準偏差であり、ΔＡ（ｎ）はΔＧ（ｎ）の絶対値平均である。９６個の細胞集合体１４の代表として、１個の細胞集合体１４の結果を、図１３に示す。

図１３は、試験例２に係るΔＳ（ｎ）とΔＡ（ｎ）の時間変化を示すグラフである。図１３において、実線がΔＳ（ｎ）の時間変化を示し、破線がΔＡ（ｎ）の時間変化を示す。図１３に示す実線と図１３に示す破線とを比較すれば明らかなように、ΔＡ（ｎ）の波形は、ΔＳ（ｎ）の波形と同じ傾向を示していた。従って、ΔＡ（ｎ）も、ΔＳ（ｎ）と同様に、細胞集合体１４の動きの速さを表すパラメータとして適切であることが分かる。

試験例３では、約２００個の細胞１５で構成される細胞集合体１４を１００個同時に撮像した動画像５０からＧ（ｎ）およびΔＧ（ｎ）を算出し、細胞集合体１４ごとにΔＳ（ｎ）を算出した。ΔＳ（ｎ）はΔＧ（ｎ）の標準偏差である。また、試験例３では、ΔＳ（ｎ）の算出に用いたのと同じ動画像５０を用いて、細胞集合体１４ごとに、特許文献２に記載のΔＰｘｙの平均値（つまり、平均ピクセル強度）を算出した。９６個の細胞集合体１４の代表として、１個の細胞集合体１４の結果を、図１４〜図１６に示す。

図１４は、試験例３に係るΔＳ（ｎ）の時間変化を示すグラフである。図１５は、試験例３に係るΔＰｘｙの平均値の時間変化を示すグラフであって、基準フレームとして９５番目のフレームを選択した時のグラフである。９５番目のフレームを撮像した時刻は、図１４および図１５に示すＴ（９５）である。時刻Ｔ（９５）では、細胞集合体１４は弛緩完了後、収縮開始前の状態である。図１６は、試験例３に係るΔＰｘｙの平均値の時間変化を示すグラフであって、基準フレームとして１番目のフレームを選択した時のグラフである。１番目のフレームを撮像した時刻は、図１４および図１５に示すＴ（１）である。時刻Ｔ（１）では、細胞集合体１４は収縮中である。図１５および図１６から明らかなように、上記特許文献２の方法では、基準フレームとして何番目のフレームを選択するかによって、異なる波形が得られることが分かる。

以上、本開示に係る画像処理方法、画像処理装置、および動作解析システムの実施形態などについて説明したが、本開示は上記実施形態などに限定されない。特許請求の範囲に記載された範疇内において、各種の変更、修正、置換、付加、削除、および組合わせが可能である。それらについても当然に本開示の技術的範囲に属する。

１０ 動作解析システム
１４ 細胞集合体
１５ 細胞
２０ 撮像装置
２１ 撮像素子
２２ 受光素子
３０ 画像処理装置
３５ 画像取得部
３６ 解析領域設定部
３７ 差分算出部
３８ 速度算出部
４０ 画像表示装置
５０ 動画像
５１ フレーム
５２ 画素
５３ 解析領域

Claims (5)

1. 複数の細胞集合体を撮像したフレームが撮像された順番で並ぶ動画像を取得する工程と、
   前記動画像における前記細胞集合体の動作を解析する解析領域を、前記細胞集合体ごとに設定する工程と、
   前記解析領域を構成する画素ごとにｎ（ｎは１以上の自然数）番目のフレームにおける階調値とｎ＋ｍ（ｍは１以上の自然数）番目のフレームにおける階調値との差分を算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程と、
   前記細胞集合体の動く速さを表すパラメータとして前記差分の標準偏差および前記差分の絶対値平均から選ばれる少なくとも１つを前記細胞集合体ごとに算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う工程とを有する、画像処理方法。
2. 前記細胞集合体の動く速さを表すパラメータとして前記差分の標準偏差および前記差分の絶対値平均から選ばれる少なくとも１つを前記細胞集合体ごとに算出することを、ｍを２以上の自然数に固定すると共にｎを変更しながら繰り返し行うことにより、ｎが１つずつ異なる算出値の集合を得る工程を有する、請求項１に記載の画像処理方法。
3. 複数の細胞集合体を撮像したフレームが撮像された順番で並ぶ動画像を取得する画像取得部と、
   前記動画像における前記細胞集合体の動作を解析する解析領域を、前記細胞集合体ごとに設定する解析領域設定部と、
   前記解析領域を構成する画素ごとにｎ（ｎは１以上の自然数）番目のフレームにおける階調値とｎ＋ｍ（ｍは１以上の自然数）番目のフレームにおける階調値との差分を算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う差分算出部と、
   前記細胞集合体の動く速さを表すパラメータとして前記差分の標準偏差および前記差分の絶対値平均から選ばれる少なくとも１つを前記細胞集合体ごとに算出することを、ｍを固定すると共にｎを変更しながら、繰り返し行う速度算出部とを有する、画像処理装置。
4. 前記速度算出部は、前記細胞集合体の動く速さを表すパラメータとして前記差分の標準偏差および前記差分の絶対値平均から選ばれる少なくとも１つを前記細胞集合体ごとに算出することを、ｍを２以上の自然数に固定すると共にｎを変更しながら繰り返し行うことにより、ｎが１つずつ異なる算出値の集合を得る、請求項３に記載の画像処理装置。
5. 請求項３または４に記載の画像処理装置と、
   前記画像処理装置で処理される前記動画像を撮像する撮像装置とを備える、動作解析システム。