JP6402717B2 - 画像解析装置、画像解析方法、画像解析プログラム、細胞の製造方法、細胞の製造装置、細胞の培養方法、および細胞の培養装置 - Google Patents

Description

本発明は、画像解析装置、画像解析方法、画像解析プログラム、細胞の製造方法、細胞の製造装置、細胞の培養方法、および細胞の培養装置に関する。
本願は、２０１３年９月１８日に出願された日本国特願２０１３−１９３５４５号及び２０１３年９月１８日に出願された日本国特願２０１３−１９３５２０号に基づき優先権を主張し、その内容をここに援用する。

近年、幹細胞を再生医療や医薬品開発へ利用する試みが行われている。幹細胞は、多分化能と増殖能を兼ね備えているため、理論上は、任意の細胞を任意の量、製造することができる。この技術の応用例として、心筋細胞の培養がある。心筋細胞は、終末分化細胞であるため、虚血性心疾患などにより傷害を受けると再生が難しい。そのため、幹細胞を心筋細胞に分化させて、得られた心筋細胞を心臓に移植する治療法が考えられている。

しかしながら、幹細胞に対して分化誘導を行ったとしても、目的の細胞に分化する細胞は、必ずしも多くない。そのため、分化後の細胞を利用したい場合には、分化した細胞と分化しなかった細胞を識別し、分離する必要がある。実際の応用を考えた場合には、安全性の確保や製造効率の面から、この細胞の選別の過程は、自動で行われることが望ましい。また、分化した細胞を傷害しないよう、非侵襲的な手法で行われることが望ましい。

時系列画像の解析技術の一つとして、画像内の物体の動きを追跡する技術がある。この技術を応用して、心臓や心筋細胞を時系列で観察した画像から心臓や心筋細胞の拍動を測定する技術が、非特許文献１及び非特許文献２に開示されている。培養心筋細胞は、自律的な拍動を行うという特性を持っている。そこで、これらの技術と培養心筋細胞の特性を利用することで、自動的かつ非侵襲的に、幹細胞から分化した心筋細胞を識別することができる。

S.Loncaric, Z.Majcenic, "Optical flow algorithm for cardiac motion estimation", Annual Conference of the IEEE Engineering in Medicine and Biology Society, 2000, 22nd, Vol.1, pp415-417 A.Kamgoue, J.Ohayon, Y.Usson, L.Riou, P.Tracqui, "Quantification of cardiomyocyte contraction based on image correlation analysis", Cytometry. Part A : The Journal of the International Society for Advancement of Cytometry, 2009, Vol.75, pp298-308

しかしながら、従来の技術では、心筋細胞の拍動だけでなく、観察画像内の他の画像要素の動きを心筋細胞の拍動として認識してしまう場合があり、心筋細胞の拍動のような自律的な運動を行う画像要素と、それ以外の画像要素とを、正確に識別することができない場合があった。

また、従来の技術では、一つの心筋細胞による動きが複数の小領域に分離されて検出されたり、逆に、複数の心筋細胞が1つの領域として検出されてしまう場合があった。すなわち、同一の心筋細胞（運動物体）を撮像した画像要素を正確にグループ化することができない場合があった。

本発明の一態様は、時系列画像から、自律的な運動を行う画像要素を、より正確に識別することを目的とする。

また、本発明の別の態様は、同一の運動物体を撮像した可能性がある一群の画像要素を、より正確にグループ化することを目的とする。

本発明の一態様の画像解析装置は、拍動する心筋細胞の時系列画像から画像要素の動きの大きさを算出する第１の算出部と、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出部と、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出部と、前記第２の算出部により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別部と、を備える画像解析装置である。

また、本発明の一態様の画像解析方法は、コンピュータが、拍動する心筋細胞の時系列画像から画像要素の動きの大きさを算出し、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出し、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出し、前記算出された特徴量に基づいて前記心筋細胞の拍動領域を識別する、画像解析方法である。

また、本発明の一態様の画像解析プログラムは、コンピュータに、拍動する心筋細胞の時系列画像から画像要素の動きの大きさを算出させ、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出させ、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出させ、前記算出された特徴量に基づいて前記心筋細胞の拍動領域を識別させる、画像解析プログラムである。

また、本発明の一態様の細胞の製造方法は、拍動する心筋細胞の時系列画像を撮像する撮像手順と、前記時系列画像から画像要素の動きの大きさを算出する第１の算出手順と、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出手順と、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出手順と、前記第２の算出手順により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別手順と、を有する細胞の製造方法である。

また、本発明の一態様の細胞の製造装置は、拍動する心筋細胞の時系列画像を撮像する撮像部と、前記時系列画像から画像要素の動きの大きさを算出する第１の算出部と、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出部と、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出部と、前記第２の算出部により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別部と、を有する細胞の製造装置である。

また、本発明の一態様の細胞の培養方法は、拍動する心筋細胞の時系列画像を撮像する撮像手順と、前記時系列画像から画像要素の動きの大きさを算出する第１の算出手順と、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出手順と、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出手順と、前記第２の算出手順により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別手順と、を有する細胞の培養方法である。

また、本発明の一態様の細胞の培養装置は、拍動する心筋細胞の時系列画像を撮像する撮像部と、前記時系列画像から画像要素の動きの大きさを算出する第１の算出部と、前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出部と、前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出部と、前記第２の算出部により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別部と、を有する細胞の培養装置である。

本発明の一態様によれば、時系列画像から、自律的な運動を行う画像要素を、より正確に識別することができる。

また、本発明の別の態様によれば、同一の運動物体を撮像した可能性がある一群の画像要素を、より正確にグループ化することができる。

画像解析部（画像解析装置）を含むインキュベータを正面から見た透視図の一例である。 インキュベータを上面から見た透視図の一例である。 第１実施形態に係る画像解析部の機能構成の一例を示す図である。 細胞の観察、撮像、及び画像処理の流れを示すフローチャートの一例である。 第１実施形態に係る領域抽出部により実行される、時系列画像からの運動領域の抽出処理の流れを示すフローチャートの一例である。 撮像部により撮像された画像の一部に含まれる細胞と、画像から領域抽出部によって抽出された運動領域との一例を示す図である。 領域抽出部によって抽出された各運動領域の動きの大きさの推移を示すグラフの一例である。 各運動領域の動きの大きさの変化率の推移を示すグラフの一例である。 各運動領域の動きの大きさの自己相関の推移を示すグラフの一例である。 各運動領域の動きの大きさの周波数スペクトルの一例を示すグラフの一例である。 識別部による運動領域の種類の識別手法を説明するための図である。 第２実施形態に係る画像解析部の機能構成の一例を示す図である。 撮像部により撮像された画像の一部に含まれる心筋細胞と、画素抽出部により「動きが大きい」と判定された画素、およびグループ化部によって抽出された運動領域の一例を示す図である。 第３実施形態に係る画像解析部の機能構成の一例を示す図である。 第４実施形態に係る画像解析部の機能構成の一例を示す図である。 領域再設定部による処理の流れの一例を模式的に示す図である。

以下、図面を参照して、本発明の画像解析装置、画像解析方法、画像解析プログラム、細胞の製造方法、および細胞の製造装置の実施の形態について説明する。

＜第１実施形態＞
［インキュベータの構成］
図１は、画像解析部（画像解析装置）を含むインキュベータ１０を正面（図中、＋Ｙ方向）から見た透視図の一例である。また図２は、インキュベータ１０を上面（図中、＋Ｚ方向）から見た透視図の一例である。図１および図２においては、互いに直交するＸＹＺ座標系を用いて位置関係を示している。

インキュベータ１０は、上部ケーシング１２と下部ケーシング１３とを有している。上部ケーシング１２は、下部ケーシング１３の上に載置される。なお、上部ケーシング１２と下部ケーシング１３との内部空間は、ベースプレート１４によって上下に仕切られている。

まず、上部ケーシング１２の構成の概要について説明する。上部ケーシング１２の内部には、細胞の培養を行う恒温室１５が形成されている。この恒温室１５は、図示しない温度調整部及び湿度調整部を有しており、恒温室１５内は、細胞の培養に適した環境（例えば温度３７℃、湿度９０％の雰囲気）に維持されている。

恒温室１５の前面には、大扉１６、中扉１７、小扉１８が配置されている。大扉１６は、上部ケーシング１２及び下部ケーシング１３の前面を覆っている。中扉１７は、上部ケーシング１２の前面を覆っており、大扉１６の開放時に恒温室１５と外部との環境を隔離する。小扉１８は、細胞を培養する培養容器１９を搬出入するための扉であり、中扉１７に取り付けられている。この小扉１８から培養容器１９を搬出入することで、恒温室１５の環境変化を抑制することが可能となる。なお、大扉１６、中扉１７、小扉１８は、パッキンＰ１、Ｐ２、Ｐ３によりそれぞれ気密性が維持されている。また、恒温室１５には、ストッカー２０、容器搬送部３０、観察ユニット４０、載置台４５が配置されている。載置台４５は、小扉１８の手前に配置されており、培養容器１９は、小扉１８を介して搬出入される。

ストッカー２０は、上部ケーシング１２の前面（＋Ｙ方向）からみて恒温室１５の左側に配置される。ストッカー２０は、複数の棚を有しており、ストッカー２０の各々の棚には培養容器１９を複数収納することができる。各々の培養容器１９は、培養の対象となる幹細胞を培地とともに収容している。また、幹細胞には心筋細胞への分化誘導が行われており、一部の細胞が心筋細胞に分化している。心筋細胞に分化した細胞は、自律的な拍動を行っている。

容器搬送部３０は、上部ケーシング１２の前面からみて恒温室１５の中央に配置される。この容器搬送部３０は、ストッカー２０から観察ユニット４０の試料台４１に培養容器１９を運搬したり、載置台４５からストッカー２０や観察ユニット４０に培養容器１９を運搬したり、観察ユニット４０の試料台４１からストッカー２０に培養容器１９を運搬したりする。

容器搬送部３０は、多関節アームを有する垂直ロボット３１と、回転ステージ３２と、ミニステージ３３と、アーム部３４とを有している。回転ステージ３２は、垂直ロボット３１の先端部に回転軸３５を中心として水平方向（Ｘ−Ｙ平面方向）に１８０°回転可能に取り付けられている。そのため、回転ステージ３２は、ストッカー２０、試料台４１及び載置台４５に対して、アーム部３４をそれぞれ対向させることができる。ミニステージ３３は、回転ステージ３２に対して水平方向に摺動可能に取り付けられている。ミニステージ３３には、培養容器１９を把持するアーム部３４が取り付けられている。

観察ユニット４０は、上部ケーシング１２の前面からみて恒温室１５の右側に配置される。この観察ユニット４０は、培養容器１９内の細胞の画像を一定時間間隔毎に撮像して観察することができる。以下、このような観察を「タイムラプス観察」とも称する。

観察ユニット４０は、上部ケーシング１２のベースプレート１４の開口部に嵌め込まれて配置される。観察ユニット４０は、試料台４１と、試料台４１の上方に張り出したスタンドアーム４２と、位相差観察用の顕微光学系及び撮像部１６０（後述）を内蔵した本体部４３とを有している。試料台４１及びスタンドアーム４２は恒温室１５に配置される一方で、本体部４３は下部ケーシング１３内に収納される。

試料台４１は、透光性の材質で構成されており、その上に培養容器１９を載置することができる。この試料台４１は、水平方向に移動可能に構成されており、上面に載置した培養容器１９の位置を調整することができる。また、スタンドアーム４２にはＬＥＤ光源が内蔵されている。そして、撮像部１６０は、スタンドアーム４２によって試料台４１の上側から透過照明された培養容器１９の細胞を、顕微光学系を介して撮像することで細胞の顕微鏡画像を取得することができる。

次に、下部ケーシング１３の構成の概要を説明する。下部ケーシング１３の内部には、本体部４３とインキュベータ１０の制御部５０が収納されている。制御部５０は、ＣＰＵ（Central Processing Unit）等のプロセッサ、ＨＤＤ（Hard Disk Drive）、フラッシュメモリ、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）等の記憶部１５０（後述）などを備えている。記憶部１５０には、ストッカー２０に収納されている各培養容器１９に関する管理データや、撮像部１６０により撮像された顕微鏡画像のデータなどが記憶されている。また、記憶部１５０には、画像処理の処理結果や処理の途中で生じる一時的なデータや、ＣＰＵによって実行されるプログラムが記憶される。

なお、上記の管理データには、個々の培養容器１９の識別情報であるインデックスデータ、ストッカー２０における培養容器１９の収納位置、培養容器１９の種類及び形状（ウェルプレート、ディッシュ、フラスコなど）、培養容器１９で培養されている細胞の種類、培養容器１９の観察スケジュール、タイムラプス観察時の撮像条件（対物レンズの倍率、容器内の観察地点等）、などが含まれている。また、ウェルプレートのように複数の小容器で同時に細胞を培養できる培養容器１９については、各々の小容器ごとにそれぞれ管理データが生成されてよい。

制御部５０は、温度調整部、湿度調整部、容器搬送部３０及び観察ユニット４０と接続されている。制御部５０は、所定のプログラムに従ってインキュベータ１０の各部を統括的に制御する。一例として、制御部５０は、温度調整部及び湿度調整部をそれぞれ制御して恒温室１５内を所定の環境条件に維持する。また、制御部５０は、所定の観察スケジュールに基づいて、容器搬送部３０、および観察ユニット４０を制御して、培養容器１９内の細胞のタイムラプス観察を実行する。さらに、制御部５０の画像解析部１００は、タイムラプス観察で取得された時系列画像を解析し、自律的な拍動を行う心筋細胞を識別する。

［画像解析部の構成］
以下、自律的な拍動を行う心筋細胞の識別について説明する。図３は、第１実施形態に係る画像解析部１００の機能構成の一例を示す図である。画像解析部１００は、例えば、第１の算出部１１０と、領域抽出部１２０と、第２の算出部１３０と、識別部１４０と、記憶部１５０とを備える。これらの機能部は、例えば、制御部５０のＣＰＵが記憶部１５０に格納されたプログラムを実行することにより機能するソフトウェア機能部である。また、これらの機能部の一部または全部は、ＬＳＩ（Large Scale Integration）やＡＳＩＣ（Application Specific Integrated Circuit）等のハードウェア機能部であってもよい。また、図１では、画像解析部１００が制御部５０の一機能であるかのように示したが、画像解析部１００（画像解析装置）は、インキュベータ１０から独立したコンピュータにより実現されてもよい。

撮像部１６０は、培養されている心筋細胞を、所定周期で繰り返し撮像する。ここで、観察ユニット４０の顕微光学系の対物レンズは２倍であり、開口数は０．０９である。また、撮像部１６０は、１０００×１０００ピクセルの解像度で撮像を行う。顕微光学系の対物レンズが２倍であるため、撮像部１６０は、１回の撮像で広い範囲の細胞を撮像することができる。そのため、画像解析部１００は、心筋細胞に分化し、拍動を行っている細胞がわずかだったとしても、これを効率よく見つけることができる。なお、上記の観察条件は一例であり、本実施形態は上記のものに限られず適用することができる。撮像部１６０は、撮像した時系列画像を第１の算出部１１０に出力する。

第１の算出部１１０は、撮像部１６０から入力された時系列画像のうち、例えば、時間的に隣接する（１つおき、２つおき等でもよい）画像の組を複数取り出し、各組の画像を比較して、画像の画素毎に動きの大きさ（後述）を算出する。第１の算出部１１０は、算出した各時系列画像についての画素毎の動きの大きさを、領域抽出部１２０に出力する。

領域抽出部１２０は、第１の算出部１１０から入力された各時系列画像についての画素毎の動きの大きさに基づいて、画像内の動きのある領域（以下、「運動領域」と称する）を抽出する。領域抽出部１２０は、各運動領域を示す識別情報と、運動領域の動きの大きさとの組を、第２の算出部１３０に出力する。

ここで、動きの大きさとは、例えば、ある画像（画像Ａ）と、その直後の画像（画像Ｂ）との間で、物体の同じ箇所を映したと考えられる画素の対応付けを行い、画像Ｂの画素毎に、画像Ａにおける対応する画素との距離を求めたものである。対応付けは、例えば、公知のオプティカルフローやブロックマッチング法等の手法を用いて行うことができる。
動きの大きさは、時系列画像のうち時間的に隣接する画素の組について得られるため、時系列画像のフレーム数よりも１少ない数だけ算出される。

第２の算出部１３０は、領域抽出部１２０から入力された各運動領域の動きの大きさから、例えば複数の特徴量を算出する。第２の算出部１３０は、算出した特徴量を識別部１４０に出力する。識別部１４０は、第２の算出部１３０から入力された特徴量に基づいて、各運動領域の運動が心筋細胞の拍動であるかどうかを識別する。識別部１４０は、識別結果を運動領域の識別情報と対応付ける。

［処理フロー］
図４は、細胞の観察、撮像、及び画像処理の流れを示すフローチャートの一例である。
（ステップＳ４００）まず、制御部５０が、管理データ内の観察スケジュールを参照し、観察時期が到来した培養容器１９の搬送命令を容器搬送部３０に出力する。容器搬送部３０は、制御部５０から入力された培養容器１９の搬送命令に基づき、前記培養容器１９を試料台に搬送する。

（ステップＳ４０２）次に、撮像部１６０が、培養容器１９内の細胞を所定周期で複数回撮像する。撮像部１６０は、例えば、画像を１５ｆｐｓ（フレーム／秒）で４秒間取得する。この場合、一回の観察で得られる時系列画像は、６０フレームとなる。撮像部１６０は、撮像した細胞の時系列画像を第１の算出部１１０に出力する。

（ステップＳ４０４）次に、第１の算出部１１０が、前述したように、撮像部１６０から入力された細胞の時系列画像における各画素と、その直後の画像の各画素との輝度の比較に基づいて、画素の動きの大きさを算出する。第１の算出部１１０は、算出した動きの大きさを領域抽出部１２０に出力する。

（ステップＳ４０６）次に、領域抽出部１２０が、第１の算出部１１０から入力された動きの大きさに基づいて、画像内の運動領域を抽出する。運動領域とは、動きの大きさが基準よりも大きい画素の集合である（詳しくは、後述する）。領域抽出部１２０は、抽出した運動領域のそれぞれに対して個別の識別情報を割り当て、各運動領域に含まれる画素の座標群と、識別情報とを対応付ける。

（ステップＳ４０８）次に、領域抽出部１２０は、各運動領域に属する画素の動きの大きさから、各運動領域の動きの大きさを算出する。領域抽出部１２０は、算出した各運動領域の動きの大きさと、運動領域の識別情報との組を、第２の算出部１３０に出力する。

（ステップＳ４１０）次に、第２の算出部１３０が、各運動領域の動きの大きさに基づいて、例えば、６つの特徴量を算出する。６つの特徴量は、例えば、動きの大きさの最大値と、突出度と、時間周期性と、周波数スペクトルの分散と、高周波成分の低周波成分に対する比と、中周波成分の低周波成分に対する比である。なお、これらの一部が省略されてもよいし、他の特徴量が算出および利用されてもよい。第２の算出部１３０は、算出した６つの特徴量と運動領域の識別情報との組とを識別部１４０に出力する。

（ステップＳ４１２）そして、識別部１４０が、第２の算出部１３０から入力された６つの特徴量に基づいて、各々の運動領域が心筋細胞の拍動を示しているか否かを識別する。

ここで、領域抽出部１２０による時系列画像からの運動領域の抽出処理（図４におけるステップＳ４０６の処理）について説明する。図５は、第１実施形態に係る領域抽出部１２０により実行される、時系列画像からの運動領域の抽出処理の流れを示すフローチャートの一例である。
（ステップＳ５００）まず、領域抽出部１２０は、時系列画像の各画像について、画素毎の動きの大きさを、それぞれ閾値と比較し、閾値以上のものと閾値未満のものとに二値化することで、動きの大きい画素と小さい画素とに分ける。

（ステップＳ５０２）次に、領域抽出部１２０は、時系列画像のうち時間的に連続する５フレームの画像を取り出す。領域抽出部１２０は、６０フレーム分の時系列画像＃０〜＃５９のうち、動きの大きさが求められた画像を画像＃１〜＃５９とすると、例えば、１回目は画像＃１〜＃５、２回目は画像＃２〜＃６、３回目は画像＃３〜＃７というように、順に５フレーム分の画像を取り出す（これと逆に、降順であっても構わない）。なお、時系列画像全体から５フレーム分の画像を取り出すのは、あくまで一例であり、領域抽出部１２０は、任意の数の画像を取り出してよい。

（ステップＳ５０４）次に、領域抽出部１２０は、取り出した５フレームの画像のうち少なくとも２フレームで動きが大きいと判定された画素を抽出し、記憶部１５０に記憶させる。

（ステップＳ５０６）次に、領域抽出部１２０は、画像＃１〜＃５９から全ての組み合わせの５フレームの画像を取り出したか否かを判定する。全ての組み合わせの５フレームの画像を取り出している場合（ＹＥＳ）、領域抽出部１２０は、ステップＳ５０８に遷移する。全ての組み合わせの５フレームの画像を取り出していない場合（ＮＯ）、領域抽出部１２０は、ステップＳ５００に戻り、次の５フレームの画像を取り出す。

（ステップＳ５０８）全ての組み合わせの５フレームの画像を取り出した場合、領域抽出部１２０は、各ルーチンのステップＳ５０２の処理において記憶部１５０に記憶された画素を全て参照し、記憶部１５０に記憶された画素が互いに隣接する領域を、一つの運動領域とする。領域抽出部１２０は、各運動領域に個別の識別情報を割り当て、前記運動領域が占める座標と合わせて記憶部１５０に記憶させる。

図６は、撮像部１６０により撮像された画像の一部（６００）に含まれる細胞６１０（１）〜６１０（３）と、画像から領域抽出部１２０によって抽出された運動領域ａ〜ｃとの一例を示す図である。運動領域としては、例えば、以下に説明する３種類のものが抽出される。運動領域ａは、心筋細胞の拍動に基づく運動が抽出された運動領域である。以下、係る運動領域を「拍動領域」と称する。運動領域ｂは、培養液の液揺れによって抽出される運動領域である。以下、係る運動領域を「液揺れ領域」と称する。運動領域ｃは、例えば、その他のノイズ要素によって抽出される運動領域である。以下、係る運動領域を「ノイズ領域」と称する。

ここで、運動領域毎の動きの大きさの算出方法（図４におけるステップＳ４０８）について詳述する。領域抽出部１２０は、以下のように、各運動領域内の画素毎の動きの大きさを、時系列画像全てについて確認する。まず、領域抽出部１２０は、運動領域内の画素のうち最も大きな動きの大きさを示した画素（「領域内最大運動画素」）を抽出する。次に、領域抽出部１２０は、領域内最大運動画素の周囲５×５ピクセルの画像領域を選択する。そして、領域抽出部１２０は、選択した５×５ピクセルの画像領域の動きの大きさの平均値を、時系列画像すべてについて算出する。領域抽出部１２０は、算出された運動領域毎の動きの大きさを、対応する運動領域の識別情報と対応付けて第２の算出部１３０に出力する。

［各運動領域の動きの大きさの特徴の概要］
ここで、各運動領域の動きの大きさの特徴について説明する。図７は、領域抽出部１２０によって抽出された各運動領域の動きの大きさの推移を示すグラフの一例である。図７のグラフの横軸は時系列画像のフレームであり、縦軸は動きの大きさの値である。また、図７において、（ａ）は拍動領域を、（ｂ）は液揺れ領域を、（ｃ）はノイズ領域を示す。

（ａ）心筋細胞の拍動は、自律的な収縮（運動）状態と弛緩（静止）状態の繰り返しである。従って、図７（ａ）に示すように、拍動領域の動きの大きさは、周期性のあるものとなる。本実施形態のサンプリング条件において、心筋細胞が収縮を始めてから弛緩するまでの時間の分布は、概ね１回当り５〜１０フレーム分に収まる。一方、心筋細胞が静止している時間は、細胞によって大きく異なっている。心筋細胞には、静止している時間がほとんどなく、短い周期で収縮を繰り返す細胞もあれば、例えば４秒間の観察で一度しか収縮しない細胞もある。

（ｂ）培養液の液揺れは、容器搬送部３０が培養容器１９を試料台４１上に輸送するときの振動により、培養容器１９内の培養液が揺り返されることで生じる。培養液の揺り返しの周期は、ほとんどの場合、心筋細胞の拍動の周期と比較して短い。従って、培養液の揺り返しによって動かされる浮遊する細胞やゴミなどの動きも心筋細胞の拍動の周期と比較して概ね短い。図７（ｂ）に示すように、液揺れ領域の動きの大きさは、拍動領域よりも短い周期で変化する。

（ｃ）ノイズ領域は、あるフレームのある画素で発生したノイズと別のフレームの前記画素の近傍で発生したノイズとが、同じ画像要素が動いたものと誤って判定されることで抽出されたものである。実際には、同じ画像要素が動いたわけではないため、検出される動きの大きさは、時間的な連続性に乏しい。従って、図７（ｃ）に示すように、ノイズ領域の動きの大きさには、周期性があまり見られないことが多い。

以上のように、運動領域の動きの大きさは、運動領域の種類に応じてそれぞれ特徴を持つ。しかしながら、実際に算出される運動領域の動きの大きさは、個々の運動領域によって差が大きく、単一の尺度で運動領域の種類を識別することは難しい。そこで、本実施形態の画像解析部１００では、動きの大きさに基づく複数の特徴量を算出し、これらを用いて総合的に運動領域の種類を識別する。

［特徴量の算出］
特徴量の算出について詳述する。本実施形態では、第２の算出部１３０は、運動領域の動きの大きさに基づいて、以下の６つの特徴量を算出する。第２の算出部１３０は、算出した６つの特徴量と領域の識別情報とを対応付けて、識別部１４０に出力する。

（１）「動きの大きさの最大値」は、運動領域毎の動きの大きさの最大の値である。第２の算出部１３０は、以下の式（１）により「動きの大きさの最大値」を算出する。式中、「ｔ」は時間、すなわち時系列画像のフレームの値を、「Ｄ（ｔ）」は動きの大きさの値を表す。例えば、ある運動領域の動きの大きさを、各フレームごとにＤ（１）〜Ｄ（５９）とすると、「動きの大きさの最大値」はＤ（１）〜Ｄ（５９）のうちの最も大きな値である。動きの大きさの最大値を特徴量として求めることによって、心筋細胞ではなく、極端に動きが少ない運動領域などを適切に除外することができる。

（２）「突出度」は、基準フレームを変えながら算出される、動きの大きさの変化の割合（以下、「変化率」と称する）の最大値である。変化率は、各運動領域について、着目した基準フレームの前の所定数のフレームにおける動きの大きさの総和と、基準フレームの後の所定数のフレームにおける動きの大きさの総和との比較に基づき算出される値である。従って、基準フレームの前の所定数のフレームで動きの大きさの値が小さく、基準フレームの後の所定数のフレームで動きの大きさの値が大きいとき、変化率は大きな値になる。一方、基準フレームの後の所定数のフレームで動きの大きさが大きくても、基準フレームの前の所定数のフレームで動きの大きさが大きければ、変化率は大きくならない。以下、上記「所定数」を、窓サイズと称する。第２の算出部１３０は、例えば、式（２）により「突出度」を算出する。式（２）中、ｗは窓サイズである。

ある運動領域の動きの大きさＤ（１）〜Ｄ（５９）について突出度を算出するとき、第２の算出部１３０は、窓サイズｗが５であれば、まず第６フレームを基準フレームとする。第２の算出部１３０は、基準フレームの５フレーム前である第１フレームから基準フレームまでの動きの大きさＤ（１）〜Ｄ（６）の総和を算出する。また、第２の算出部１３０は、基準フレームから基準フレームの５フレーム後である第１１フレームまでの動きの大きさＤ（６）〜Ｄ（１１）の総和を算出する。

第２の算出部１３０は、算出した動きの大きさＤ（６）〜Ｄ（１１）の総和と、動きの大きさＤ（１）〜Ｄ（６）の総和との比較に基づき（例えば、Ｄ（６）〜Ｄ（１１）の総和をＤ（１）〜Ｄ（６）の総和で除算することで）、変化率を算出する。そして、第２の算出部１３０は、例えば、１回目は第６フレーム、２回目は第７フレームというように、順に基準フレームを変えながら、第５４フレームを基準とした場合まで、変化率を算出する。突出度は、各運動領域について、基準フレームを変えながら算出された変化率のうちの最大の値である。

図８は、各運動領域の動きの大きさの変化率の推移を示すグラフの一例である。前記グラフの横軸は基準フレームであり、前記グラフの縦軸は変化率である。また、図８において、（ａ）は拍動領域を、（ｂ）は液揺れ領域を、（ｃ）はノイズ領域を示す。

心筋細胞の拍動には、弛緩（静止状態）と収縮（運動状態）が含まれるため、静止状態から運動状態に移るとき、動きの大きさの変化率は大きくなる。従って、図８（ａ）に示すように、拍動領域の動きの大きさの変化率のグラフは、特に静止状態がある程度長く（例えば５フレーム程度以上あり）、その後、大きく立ち上がったような場合に、明確なピークを示す。

一方、液揺れ領域のように静止状態がほとんどない運動領域では、動きの大きさの変化率は大きくならない。従って、図８（ｂ）に示すように、液揺れ領域の動きの大きさの変化率のグラフは、明確なピークを示さない。

また、ノイズ領域では、静止状態はあっても、窓サイズの５フレームという時間で持続して動くことはあまりないため、動きの大きさの変化率は小さくなる。従って、図８（ｃ）に示すように、ノイズ領域の動きの大きさの変化率のグラフは、明確なピークを示さない。

このように、突出度は、静止状態が明確に現れる拍動領域の識別に有効な特徴量である。また、突出度は、比較的長時間（例えば４秒など）に一度しか収縮せず、本実施形態における観察条件では周期性を確認できない心筋細胞を識別するためにも、有効な特徴量である。

ここで、心筋細胞が１回の拍動に要する時間の分布が、概ね５〜１０フレーム程度に収まるとすると、窓サイズを３〜４フレーム以上とすることで、心筋細胞の拍動を効果的に検出することができる。「概ね収まる」とは、例えば、９０［％］程度が収まること、正規分布にして−２σ〜＋１σ程度が収まること、その他の基準で定義される。一方、窓サイズをある程度大きく設定することで、短い周期で周期的な運動を行っている領域の突出度を低く評価するフィルターの効果を強めることができる。フィルター効果を高めるためには、例えば、窓サイズを１０フレーム程度としてもよい。検出感度とフィルター効果を両立させるためには、窓サイズを、心筋細胞が１回の拍動に要する時間の分布が概ね収まる範囲（５〜１０フレーム）の下限値の１／２〜上限値程度とすると好適である。本実施形態では、その一例として５フレームを窓サイズとしている。これによって、画像解析部１００は、心筋細胞の拍動を、より正確に、且つ漏れなく抽出することができる。

なお、突出度の算出式は、式（２）に限定されない。例えば、式（２）において、除算しているところを減算に置き換え、その絶対値を算出してもよい。以上のように、突出度の算出は、特定の式や方法に限定されず、動きの大きさの変化の程度を反映するものであればよい。

（３）「時間周期性」は、例えば、窓間隔を変えながら算出される、動きの大きさの自己相関の最大値である。第２の算出部１３０は、以下の式（３）により「時間周期性」を算出する。式中、τは、窓間隔である。また、Ｎは、総フレーム数（例えば５９）である。第２の算出部１３０は、各運動領域についてフレーム毎に自己相関を求め、その最大値を時間周期性として算出する。

第２の算出部１３０による「時間周期性」の算出について、より詳細に説明する。例えば、窓間隔τが１であるとき、第２の算出部１３０は、動きの大きさＤ（１）〜Ｄ（５８）と、１フレームずれた動きの大きさＤ（２）〜Ｄ（５９）とを、１対１で対応付ける。
第２の算出部１３０は、対応付けられた動きの大きさの組についてそれぞれ積（Ｄ（１）×Ｄ（２）、Ｄ（２）×Ｄ（３）、‥）を算出する。次に、第２の算出部１３０は、算出された積の平均値を算出する。算出される積は、組として対応付けられた動きの大きさが共に大きいときに大きい値となる。また、平均値であるので、動きの大きさのフレーム毎の変化について、元の動きの大きさと窓間隔分ずらした動きの大きさとが組全体において類似しているときに、自己相関は大きくなる。そして、第２の算出部１３０は、例えば、１回目は１フレーム、２回目は２フレームというように順に窓間隔を変えていき、５８フレームまでフレームをずらしたときの自己相関を算出する。「時間周期性」は、各運動領域について、フレーム毎に算出された自己相関の最大値である。従って、動きの大きさが周期性を持つとき、「時間周期性」は大きな値を示す。その逆に、動きの大きさが周期性に乏しいとき、「時間周期性」は低い値を示す。

図９は、各運動領域の動きの大きさの自己相関の推移を示すグラフの一例である。前記グラフの横軸は時系列画像の窓間隔τであり、前記グラフの縦軸は自己相関である。また、図９において、（ａ）は拍動領域を、（ｂ）は液揺れ領域を、（ｃ）はノイズ領域を示す。前述したように、拍動領域や液揺れ領域の運動は、ある程度の周期性を示す。そのため、図９（ａ）及び（ｂ）に示すように、拍動領域や液揺れ領域の自己相関は窓間隔の変化に応じたピークを示す。一方、ノイズ領域の運動は周期的ではないため、図９（ｃ）に示すように、ノイズ領域の自己相関は、明確なピークを示さない。従って、「時間周期性」は、ノイズ領域を適切に除外するために有効な特徴量である。

残りの３つの特徴量（（４）〜（６））は、動きの大きさの周波数スペクトルから算出される。第２の算出部１３０は、動きの大きさの周波数スペクトルを、例えばフーリエ変換により算出する。フーリエ変換による周波数スペクトルの算出は、公知の技術であるため、ここでは説明を省略する。そして、第２の算出部１３０は、算出した動きの大きさの周波数成分を、高周波成分、中周波成分、低周波成分に分類する。本実施形態の観察条件において、「高周波成分」とは、例えば周波数１５〜２９［Ｈｚ］の周波数成分であり、「中周波成分」とは、例えば周波数７〜１２［Ｈｚ］の周波数成分であり、「低周波成分」とは、例えば周波数１〜６［Ｈｚ］の周波数成分である。図１０は、各運動領域の動きの大きさの周波数スペクトルの一例を示すグラフの一例である。前記グラフの横軸は周波数であり、前記グラフの縦軸はパワースペクトル密度である。また、図１０において、（ａ）は拍動領域を、（ｂ）は液揺れ領域を、（ｃ）はノイズ領域を示す。

（４）「高／低周波成分比」は、動きの大きさの低周波成分に対する高周波数成分の比である。図１０（ａ）及び（ｂ）に示すように、液揺れ領域の動きの大きさは、拍動領域の動きの大きさと比較して、高周波成分が多い傾向にある。従って、「高／低周波成分比」は、拍動領域と液揺れ領域とを識別するのに有効な特徴量である。
（５）「中／低周波成分比」は、動きの大きさの低周波成分に対する中周波数成分の比である。図１０（ａ）及び（ｂ）に示すように、拍動領域の動きの大きさは、液揺れ領域の動きの大きさと比較して、中周波成分が多い傾向にある。従って、「中／低周波成分比」は、拍動領域と液揺れ領域との識別に有効な特徴量である。

なお、周波数成分の分類は前述のものに限定されない。第２の算出部１３０は、周波数成分をより細かく分類してもよいし、３つではなく２つに分けてもよい。また、上記「高／低周波成分比」や「中／低周波成分比」は、周波数スペクトルに基づいた特徴量の一例である。第２の算出部１３０は、例えば、高周波成分、中周波成分のパワースペクトル密度そのものや、パワースペクトル密度を何らかの基準で正規化したもの、或いはこれらに演算を加えたものを、特徴量として算出してもよい。

（６）「周波数スペクトルの分散」とは、動きの大きさの周波数成分の分散である。前述したように、拍動領域及び液揺れ領域の動きの大きさは、周期的である。そのため、拍動領域及び液揺れ領域の周波数スペクトルの分散は小さくなる。一方、前述したように、ノイズ領域の動きの大きさは、周期的ではない。そのため、ノイズ領域の周波数スペクトルは、特定のピークを示さず、周波数スペクトルの分散は大きくなる。従って、「周波数スペクトルの分散」は、ノイズ領域を好適に除外するために有効な特徴量である。

なお、第２の算出部１３０は、前述した６つ以外の特徴量を算出してもよい。例えば、領域の輝度も、細胞の識別においては有効な指標になりうる。細胞領域は、液揺れ領域やノイズ領域に比較して、輝度が低くなることが多いからである。このように、第２の算出部１３０は、運動領域の動きの大きさや領域固有の値に基づいて、前述の特徴量以外にも様々な特徴量を算出してよい。

［運動領域の種類の識別］
以下、運動領域の識別方法について詳述する。識別部１４０は、第２の算出部１３０が算出した特徴量に基づいて、各運動領域を、拍動領域、液揺れ領域、またはノイズ領域のいずれかに分類する。運動領域の識別方法には、例えば、ｋ−近傍法を用いることができる。ｋ−近傍法による識別を行う場合、画像解析部１００は、予め記憶部１５０に記憶されている教師データを用いて識別を行う。

まず、識別部１４０は、第２の算出部１３０が分類未知の運動領域について算出した特徴量に基づいて、前記分類未知の運動領域のデータに、各特徴量を座標軸とする座標空間上の距離が近い教師データを、所望の数ｋだけ検出する。そして、識別部１４０は、ｋ個の教師データの領域の種類のうち、最も多い領域の種類を、分類未知の運動領域の種類であると判定する。識別部１４０は、運動領域の識別情報と、識別した運動領域の種類とを対応付けて記憶部１５０に記憶する。記憶部１５０に記憶された情報は、以降、教師データとして使用されてもよいし、使用されなくてもよい。識別部１４０は、ｋ−近傍法に限らず、サポートベクターマシンや、判別分析、決定木などの手法を利用して、運動領域の識別を行ってもよい。

図１１は、識別部１４０による運動領域の種類の識別手法を説明するための図である。
本実施形態において、識別部１４０は、６つの特徴量を用いて運動領域の種類の識別を行うが、図１１は、簡易的に２つの特徴量を用いた場合を示している。図１１は、横軸を特徴量１、縦軸を特徴量２としている。図１１には、それぞれの種類の運動領域の教師データがプロットされている。図１１にプロットされている「○」は拍動領域の教師データであり、「□」は液揺れ領域の教師データであり、「△」はノイズ領域の教師データである。「×」は、領域抽出部１２０によって抽出された分類未知の運動領域であり、第２の算出部１３０によって算出された特徴量に基づいて、図１１中の対応する位置にプロットされている。

ｋ−近傍法では、種類が未知のデータの識別が、近傍の教師データとの比較に基づいて行われる。ここで、ｋ＝５の場合について説明する。この場合、識別部１４０は、種類が未知の領域のデータと、その近傍の５つのデータとの比較を行う。種類が未知の領域のデータ「×」の近傍５つのデータは、図１１中の「×」を中心とする円領域７００に含まれている。「×」の近傍５つのデータのうち、３つは液揺れ領域であり、残りの２つが拍動領域である。従って、識別部１４０は、近傍５つの領域の種類の多数決から、「×」を示す分類未知の領域は液揺れ領域であると判定する。

［まとめ］
以上説明したように、第１実施形態の画像解析部１００によれば、第２の算出部１３０が、領域抽出部１２０により算出された運動領域毎の動きの大きさに基づいて特徴量を算出し、識別部１４０が、第２の算出部１３０により算出された特徴量に基づいて運動領域が拍動領域であるか否かを識別するため、拍動領域（すなわち自律的な運動をする運動領域＝画像要素）を、より正確に識別することができる。

また、第１実施形態の画像解析部１００によれば、第２の算出部１３０が、複数の特徴量を算出し、識別部１４０が、複数の特徴量に基づいて運動領域が拍動領域であるか否かを識別するため、１つの特徴量のみでは識別が難しい場合であっても、拍動領域を正確に識別することができる。

また、第１実施形態の画像解析部１００によれば、第２の算出部１３０が、時間的に前後する２つの窓サイズ（時間帯）間の動きの大きさの比較に基づいて突出度を算出し、識別部１４０が、突出度に基づいて運動領域が拍動領域であるか否かを識別するため、あまり高頻度に拍動しない拍動領域であっても、正確に識別することができる。また、第１実施形態の画像解析部１００によれば、第２の算出部１３０が、心筋細胞の１回の拍動に要する時間（５〜１０フレーム程度）の下限値の１／２〜上限値程度の時間を窓サイズとして突出度を算出ため、拍動領域を、更に正確に識別することができる。

また、第１実施形態の画像解析部１００によれば、第２の算出部１３０が、動きの大きさの周波数スペクトルに基づく特徴量を算出し、識別部１４０が、これを用いて運動領域が拍動領域であるか否かを識別するため、周波数特性に違いがある運動領域を正確に識別することができる。なお、拍動領域であると識別された運動領域は、例えば赤色等の色を付け、表示部に表示してもよい。また、識別された運動領域を線で囲むことで運動領域が明確になるように表示してもよい。

また、マニピュレータによって自動的に取り出された細胞のみを更にインキュベータ内で培養する培養方法も提供する。これにより、より選別された所望の細胞のみの培養を継続することができ、再生医療や創薬の研究において有効な細胞を提供することが可能となる。

さらに、特に創薬の分野において、継続された培養によって培養された所望の細胞に対して、研究対象となる薬品を添加するピペット等の添加手段を備えてもよい。添加手段を備えることで、継続して培養された細胞を用いた創薬の研究をさらに継続することも可能である。添加手段は、添加薬品のインキュベータ内への拡散によるコンタミ（雑菌等の混入）の問題を考慮し、インキュベータの外部に配置してもよい。この場合、インキュベータと添加手段を備えた装置との間を細胞が収納された培養容器を搬送する搬送装置を備えることが好ましい。

この構成によれば、薬品が添加された細胞の状態を観察し、また画像による解析も可能となる。

なお、インキュベータ１０は、拍動領域と識別された部分の細胞を自動的に取り出すマニピュレータ（取出部）を備えた、細胞の製造装置として構成されてもよい。後述する第２〜第４実施形態についても同様である。

＜第２実施形態＞
［構成］
以下、第２実施形態に係る画像解析部（画像解析装置）について説明する。第１実施形態と同様の構成要素については、説明を省略する。図１２は、第２実施形態に係る画像解析部１００Ａの機能構成の一例を示す図である。画像解析部１００Ａは、第１実施形態と同様、例えば、インキュベータ１０の制御部５０の一機能であり、或いはインキュベータ１０から独立したコンピュータにより実現される。画像解析部１００Ａは、例えば、第１の算出部１１０と、画素抽出部１２１と、グループ化部１２３と、第２の算出部１３０と、識別部１４０と、記憶部１５０とを備える。これらの機能部は、例えば、制御部５０のＣＰＵが記憶部１５０に格納されたプログラムを実行することにより機能するソフトウェア機能部である。また、これらの機能部の一部または全部は、ＬＳＩやＡＳＩＣ等のハードウェア機能部であってもよい。

第１の算出部１１０は、第１実施形態と同様に、画像の画素毎に動きの大きさを算出する。画素抽出部１２１は、第１実施形態の領域抽出部１２０と同様に、各時系列画像についての画素毎の動きの大きさに基づいて、動きが大きい画素を抽出し、記憶部１５０に記憶させる。すなわち、画素抽出部１２１は、図５のフローチャートにおけるステップＳ５００〜Ｓ５０６の処理を行う。

グループ化部１２３は、画素抽出部１２１により抽出された「動きの大きさが大きい画素」間で、動きの大きさ（動きに関する特徴量の一例である）が類似する画素を、同一の心筋細胞を撮像した一群の画素の候補としてグループ化する。グループ化部１２３によりグループ化された画素は、第１実施形態における「運動領域」と同様に扱われる。すなわち、グループ化部１２３は、各グループ（運動領域）内の画素のうち最も大きな動きの大きさを示した画素（「領域内最大運動画素」）を抽出し、領域内最大運動画素の周囲５×５ピクセルの画像領域を選択し、選択した５×５ピクセルの画像領域の動きの大きさの平均値を、時系列画像すべてについて算出する。そして、グループ化部１２３は、算出した運動領域毎の動きの大きさを、対応する運動領域の識別情報と合わせて第２の算出部１３０に出力する。

以下、グループ化処理について説明する。グループ化部１２３は、例えば、画素抽出部１２１により抽出された画素から、動きの大きさが特に大きい画素（以下、「代表画素」と称する）を幾つか抽出する。「特に大きい」とは、例えば、動きの大きさが、画素抽出部１２１が抽出する際に使用する閾値よりも大きい第２の閾値以上であることを意味する。また、これに代えて、例えば動きの大きさが上位の数［％］程度の画素を、代表画素としてもよい。そして、グループ化部１２３は、代表画素と他の全ての画素との組み合わせについて、類似度を算出し、類似度が例えば閾値以上である場合に、代表画素と他の画素を同じグループに含める。なお、グループ化部１２３は、代表画素同士を同じグループに含めてもよい。

［類似度］
類似度は、例えば、式（４）で表される正規化相互相関により算出される。式中、Ｄ（ｔ）は、動きの大きさであり、上線付きのＤは、動きの大きさの平均値であり、Ｎは、総フレーム数（例えば５９）である。これによって、動きの大きさＤ（ｔ）の推移が類似する程度を、適切に数値化することができる。

また、グループ化部１２３は、類似度の計算に画素間の距離を反映させて、距離が大きい画素の類似度が低く算出されるようにしてもよい。この場合、グループ化部１２３は、正規化相互相関等の指標値に、距離の逆数、または距離の二乗の逆数を乗算して類似度を算出してよい。式（５）は、距離の逆数を正規化相互相関に乗算して求められる類似度を示している。式中、Ｘ_１２は、画素１と画素２との間の距離である。こうすれば、距離が離れている画素であって、たまたま動きの周期が類似する画素を、同一のグループに含めてしまうのを予防することができる。

このように、類似度に基づいてグループ化された一群の画素を運動領域とすることで、より適切に運動領域を切り出すことができる。ここで、第１実施形態では、動きの大きさが大きいと判定された画素のうち、互いに隣接する画素を一群の画素、すなわち運動領域としていた。ところが、心筋拍動は、それぞれ固有の周期を持っており、同じ心筋細胞を撮像した画素であれば、同じ周期で運動するものである。また、心筋細胞の中に動きの小さい部分がある場合、その部分が、画素抽出部１２１により「動きが小さい」と判定される結果、互いに離間した画素に跨って一つの心筋細胞が撮像される場合もある。

従って、第２実施形態の画像解析部１００Ａでは、同じ周期で運動する複数の画素を、１つの運動領域として扱うために、動きの大きさを時系列の一次元信号として画素毎に抽出し、一次元信号に基づいて画素間の類似度を算出することで、同じ周期で運動する画素を抽出することができる。この結果、第２実施形態の画像解析部１００Ａは、本来一つの拍動領域に属するべき画素が別々の運動領域として扱われたり、逆に、二つの拍動領域が一つの運動領域として扱われたりする可能性を、低減することができる。図１３は、撮像部１６０により撮像された画像の一部（８００）に含まれる心筋細胞８１０（１）、８１０（２）と、画素抽出部１２１により「動きが大きい」と判定された画素ＭＰ、およびグループ化部１２３によって抽出された運動領域ｄ、ｅの一例を示す図である。図示するように、画素抽出部１２１により「動きが大きい」と判定された画素ＭＰは、同じ心筋細胞８１０を撮像したものであっても、互いに離間した画素群として現れる場合がある。第２実施形態の画像解析部１００Ａでは、このような場合でも、同じ心筋細胞を撮像した画素を、同じ運動領域に含めることができる。

グループ化部１２３により運動領域が決定されると、第２の算出部１３０および識別部１４０が、第１実施形態と同様の処理を実行する。これについての説明は省略する。

［まとめ］
以上説明した第２実施形態の画像解析部１００Ａによれば、動きの大きさの推移が類似する画素を同じグループに含めることで、同一の運動物体（心筋細胞）を撮像した可能性がある一群の画素を、より正確にグループ化することができる。

また、第２実施形態の画像解析部１００Ａによれば、相対的な距離が大きい画素の類似度が低く算出されるようにすることで、本来は同じ心筋細胞を撮像した可能性が低い画素を、同一のグループに含めてしまうのを防止することができる。

なお、第２実施形態では、グループ化部１２３は、動きが特に大きい代表画素を基準としてグループ化を行うものとしたが、これに代えて、画素抽出部１２１により抽出された画素の全ての組み合わせについて類似度を算出し、グループ化を行ってもよい。また、画素抽出部１２１が省略され、画素抽出部１２１が存在すれば動きが小さいとして除外された画素を含めた全ての画素の組み合わせについて、グループ化部１２３が類似度を算出し、グループ化を行ってもよい。

＜第３実施形態＞
以下、第３実施形態に係る画像解析部（画像解析装置）について説明する。第１実施形態と同様の構成要素については、説明を省略する。図１４は、第３実施形態に係る画像解析部１００Ｂの機能構成の一例を示す図である。画像解析部１００Ｂは、第１実施形態と同様、例えば、インキュベータ１０の制御部５０の一機能であり、或いはインキュベータ１０から独立したコンピュータにより実現される。画像解析部１００Ｂは、例えば、第１の算出部１１０と、画素抽出部１２１と、グループ化用特徴量算出部１２２と、グループ化部１２３と、第２の算出部１３０と、識別部１４０と、記憶部１５０とを備える。これらの機能部は、例えば、制御部５０のＣＰＵが記憶部１５０に格納されたプログラムを実行することにより機能するソフトウェア機能部である。また、これらの機能部の一部または全部は、ＬＳＩやＡＳＩＣ等のハードウェア機能部であってもよい。

第１の算出部１１０は、第１実施形態と同様に、画像の画素毎に動きの大きさを算出する。画素抽出部１２１は、第１実施形態の領域抽出部１２０と同様に、各時系列画像についての画素毎の動きの大きさに基づいて、動きが大きい画素を抽出し、記憶部１５０に記憶させる。すなわち、画素抽出部１２１は、図５のフローチャートにおけるステップＳ５００〜Ｓ５０６の処理を行う。

グループ化用特徴量算出部１２２は、画素抽出部１２１により「動きが大きい」と判定された各画素について、例えば複数の特徴量を算出する。グループ化用特徴量算出部１２２によって算出される特徴量は、例えば、第１実施形態で説明した、動きの大きさの変化率（式（２）の括弧内）、動きの大きさの周波数成分比等である。

そして、第２実施形態に係るグループ化部１２３は、例えば、第２実施形態と同様に、画素抽出部１２１により抽出された画素から代表画素を幾つか抽出し、代表画素と他の全ての画素との間で、複数の特徴量について、例えばｍｅａｎ ｓｈｉｆｔなどの手法によるクラスタリング分析を行い、クラスタリングによって同一のグループに属することになった画素群を、特徴量が類似する画素として一つの運動領域に設定する。ここで、特徴量には、第１の算出部１１０により算出された動きの大きさと、グループ化用特徴量算出部１２２により算出された特徴量の双方が含まれてよい。なお、グループ化部１２３は、代表画素同士を同じグループに含めてもよい。また、グループ化部１２３は、複数の特徴量について正規化相互相関を求めて類似度を算出してもよい。

このように、画素毎に算出された複数の特徴量に基づいてクラスタリング分析を行うことによりグループ化された一群の画素を運動領域とすることで、より適切に運動領域を切り出すことができる。これについては、第２実施形態と同様であるため、説明を省略する。また、グループ化部１２３により運動領域が決定されると、第２の算出部１３０および識別部１４０が、第１実施形態と同様の処理を実行する。なお、第３実施形態においても、第２実施形態と同様、相対的な距離が大きい画素が同じ運動領域に属する可能性が低くなるような処理を行ってもよい。

［まとめ］
以上説明した第３実施形態の画像解析部１００Ｂによれば、画素毎に算出された複数の特徴量に基づいてクラスタリング分析を行うことによりグループ化された一群の画素を運動領域とすることで、同一の運動物体（心筋細胞）を撮像した可能性がある一群の画素を、より正確にグループ化することができる。

なお、第３実施形態では、グループ化部１２３は、動きが特に大きい代表画素を基準としてグループ化を行うものとしたが、これに代えて、画素抽出部１２１により抽出された画素の全ての組み合わせについて類似度を算出し、グループ化を行ってもよい。また、画素抽出部１２１が省略され、画素抽出部１２１が存在すれば動きが小さいとして除外された画素を含めた全ての画素の組み合わせについて、グループ化部１２３がクラスタリング分析、または類似度の算出を行い、グループ化を行ってもよい。

＜第４実施形態＞
以下、第４実施形態に係る画像解析部（画像解析装置）について説明する。第１実施形態と同様の構成要素については、説明を省略する。図１５は、第４実施形態に係る画像解析部１００Ｃの機能構成の一例を示す図である。画像解析部１００Ｃは、第１実施形態と同様、例えば、インキュベータ１０の制御部５０の一機能であり、或いはインキュベータ１０から独立したコンピュータにより実現される。画像解析部１００Ｃは、例えば、第１の算出部１１０と、領域抽出部１２０と、第２の算出部１３０と、識別部１４０と、領域再設定部１４５と、記憶部１５０とを備える。これらの機能部は、例えば、制御部５０のＣＰＵが記憶部１５０に格納されたプログラムを実行することにより機能するソフトウェア機能部である。また、これらの機能部の一部または全部は、ＬＳＩやＡＳＩＣ等のハードウェア機能部であってもよい。第１の算出部１１０、領域抽出部１２０、第２の算出部１３０、および識別部１４０については、第１実施形態と同様であるため、説明を省略する。

領域再設定部１４５は、識別部１４０によって拍動領域であると識別された一または複数の運動領域について、同じ心筋細胞によるものであるか否かを判定し、拍動領域の再設定を行う。図１６は、領域再設定部１４５による処理の流れの一例を模式的に示す図である。図１６（Ａ）は、撮像部１６０により撮像された画像の一部（９００）に含まれる心筋細胞９１０（１）、９１０（２）と、識別部１４０によって拍動領域であると識別された拍動領域ｆ、ｇを示している。図１６（Ａ）において、同じ心筋細胞９１０（１）を撮像した筈の拍動領域ｆ、ｇは、第１の算出部１１０〜識別部１４０によって、別の拍動領域として出力されている。また、拍動領域ｇには、心筋細胞（１）とは別の心筋細胞９１０（２）が含まれてしまっている。領域再設定部１４５は、このような状況が生じ得ることに鑑み、拍動領域の再設定を行う。

まず、領域再設定部１４５は、例えば、拍動領域の形状において、くびれのある箇所を探索し、拍動領域を複数の領域に分割する。くびれのある箇所は、例えば、分割線の位置および傾きを変えながら形状を分割し、長さが極小となる分割線を求めることなどで発見することができる。また、これに代えて、領域再設定部１４５は、無作為に（ランダムに、あるいは所望のサイズの格子状に）拍動領域を分割してよい。また、領域再設定部１４５は、図１６（Ａ）に示す拍動領域ｆのように、サイズが小さい拍動領域については、分割を行わないものとしてよい。ここでは、拍動領域ｇが、拍動領域ｇ−１、ｇ−２に分割されたものとする。

次に、領域再設定部１４５は、拍動領域ｆ、ｇ−１、ｇ−２のそれぞれについて、例えば動きの大きさを算出する。領域再設定部１４５は、動きの大きさに代えて、或いは加えて、動きに関する他の特徴量を算出してもよい。拍動領域の動きの大きさは、例えば、第１実施例と同様に、拍動領域内の画素のうち最も大きな動きの大きさを示した領域内最大運動画素の周囲５×５ピクセルの画像領域を選択し、選択した５×５ピクセルの画像領域の動きの大きさの平均値として求められる。

次に、領域再設定部１４５は、拍動領域ｆ、ｇ−１、ｇ−２のそれぞれについて算出した動きの大きさ同士で、例えば相関係数を求めることにより、類似度を算出する（図１６（Ｂ）参照）。ここでは、拍動領域ｆとｇ−１の相関係数が高く、拍動領域ｇ−１とｇ−２の相関係数が低く算出されたものとする。この結果、領域再設定部１４５は、拍動領域ｆとｇ−１を、新たに一つの拍動領域ｈとして再設定し、拍動領域ｇ−２を、新たに一つの拍動領域ｉとして再設定する（図１６（Ｃ）参照）。これによって、領域再設定部１４５は、同一の心筋細胞を撮像した可能性がある一群の画素を、より正確にグループ化することができる。

［まとめ］
以上説明した第４実施形態の画像解析部１００Ｃによれば、識別部１４０によって拍動領域であると識別された一または複数の運動領域について、同じ心筋細胞によるものであるか否かを判定し、拍動領域の再設定を行うため、同一の運動物体（心筋細胞）を撮像した可能性がある一群の画素を、より正確にグループ化することができる。

以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は前述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
例えば、上記各実施形態の処理は、心筋細胞とは別の観察対象に対して行われてもよい。

また、第２実施形態または第３実施形態の機能と、第４実施形態の機能とを組み合わせてもよい。すなわち、拍動領域の識別前に、動きに関する特徴量が類似する画素（または運動領域）を処理対象の運動領域とし、拍動領域の識別後にも、動きに関する特徴量が類似する領域を拍動領域に再設定してよい。

１０‥インキュベータ、１９‥培養容器、２０‥ストッカー、３０‥容器搬送部、４０‥観察ユニット、５０‥制御部、１００、１００Ａ、１００Ｂ、１００Ｃ‥画像解析部、１１０‥第１の算出部、１２０‥領域抽出部、１２１‥画素抽出部、１２２‥グループ化用特徴量算出部、１２３‥グループ化部、１３０‥第２の算出部、１４０‥識別部、１４５‥領域再設定部、１５０‥記憶部、１６０‥撮像部

Claims (16)

1. 拍動する心筋細胞の時系列画像から画像要素の動きの大きさを算出する第１の算出部と、
   前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出部と、
   前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出部と、
   前記第２の算出部により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別部と、
   を備える画像解析装置。
2. 前記第２の算出部は、複数の特徴量を算出し、
   前記識別部は、前記複数の特徴量に基づいて前記心筋細胞の拍動領域を識別する、請求項１に記載の画像解析装置。
3. 前記第２の算出部は、前記時系列画像のうち所定の画像を基準画像として、前記画像要素の動きの大きさの変化率を算出し、前記変化率に関する情報を用いて前記特徴量を算出する、請求項１又は請求項２に記載の画像解析装置。
4. 前記第２の算出部は、前記画像要素の動きの大きさの自己相関に関する情報を用いて、前記特徴量を算出する、請求項１から請求項３のいずれか一項に記載の画像解析装置。
5. 前記第２の算出部は、前記画像要素の動きの大きさの周波数スペクトルを算出し、前記算出した周波数スペクトルに基づいて前記特徴量を算出する、請求項１から請求項４のいずれか一項に記載の画像解析装置。
6. 前記領域抽出部は、前記画像要素の動きの大きさが所定値よりも大きい画像要素を抽出し、隣接する画素要素を１つの領域とする、請求項１から請求項５のいずれか一項に記載の画像解析装置。
7. 前記領域抽出部は、前記画像要素の動きの大きさが所定値よりも大きい画像要素を抽出し、前記画像要素の動きの大きさの類似度に基づいて前記画像領域を抽出する、請求項１から請求項５のいずれか一項に記載の画像解析装置。
8. 前記領域抽出部は、前記画像要素の動きの大きさが所定値よりも大きい画像要素を抽出し、前記抽出した画像要素について算出された特徴量に基づいて前記画像領域を抽出する、請求項１から請求項５のいずれか一項に記載の画像解析装置。
9. 前記識別部が前記心筋細胞の拍動領域と識別した複数の画像領域について、同じ心筋細胞であるか否かを判定し、判定結果に基づいて、画像領域の再設定を行う領域再設定部を備える、請求項１から請求項８のいずれか一項に記載の画像解析装置。
10. 前記特徴量は、突出度、時間周期性、周波数スペクトルの分散、高周波成分の低周波成分に対する比、及び中周波成分の低周波成分に対する比のいずれかである、請求項１から請求項９のいずれか一項に記載の画像解析装置。
11. コンピュータが、
    拍動する心筋細胞の時系列画像から画像要素の動きの大きさを算出し、
    前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出し、
    前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出し、
    前記算出された特徴量に基づいて前記心筋細胞の拍動領域を識別する、画像解析方法。
12. コンピュータに、
    拍動する心筋細胞の時系列画像から画像要素の動きの大きさを算出させ、
    前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出させ、
    前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出させ、
    前記算出された特徴量に基づいて前記心筋細胞の拍動領域を識別させる、画像解析プログラム。
13. 拍動する心筋細胞の時系列画像を撮像する撮像手順と、
    前記時系列画像から画像要素の動きの大きさを算出する第１の算出手順と、
    前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出手順と、
    前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出手順と、
    前記第２の算出手順により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別手順と、
    を有する細胞の製造方法。
14. 拍動する心筋細胞の時系列画像を撮像する撮像部と、
    前記時系列画像から画像要素の動きの大きさを算出する第１の算出部と、
    前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出部と、
    前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出部と、
    前記第２の算出部により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別部と、
    を有する細胞の製造装置。
15. 拍動する心筋細胞の時系列画像を撮像する撮像手順と、
    前記時系列画像から画像要素の動きの大きさを算出する第１の算出手順と、
    前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出手順と、
    前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出手順と、
    前記第２の算出手順により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別手順と、
    を有する細胞の培養方法。
16. 拍動する心筋細胞の時系列画像を撮像する撮像部と、
    前記時系列画像から画像要素の動きの大きさを算出する第１の算出部と、
    前記画像要素の動きの大きさが所定値よりも大きい画像領域を抽出する領域抽出部と、
    前記画像領域に含まれる前記画像要素の動きの大きさに基づいて特徴量を算出する第２の算出部と、
    前記第２の算出部により算出された前記特徴量に基づいて前記心筋細胞の拍動領域を識別する識別部と、
    を有する細胞の培養装置。

Images