JP7382290B2 - 画像処理方法、プログラムおよび記録媒体 - Google Patents

Description

この発明は、細胞塊、特に内部に空洞を有する細胞塊を撮像しその構造を解析するための画像処理に関する。

例えば不妊治療を目的とした生殖補助医療においては、体外で受精させ一定期間培養した胚（受精卵）を体内に戻すことが行われる。しかしながら、その（生殖補助医療における）妊娠成功率は必ずしも高くなく、患者の精神的および経済的な負担も大きい。この問題を解決するために、培養される胚の状態を的確に判断する方法が模索されている。

例えば、胚盤胞期まで発生が進んだ胚では、内部に胞胚腔と呼ばれる空洞構造が生じ、その周囲を栄養外胚葉と呼ばれる細胞の層が覆っている。栄養外胚葉を構成する細胞の数は胚の状態を評価する指標となるものであり、したがってその数を非侵襲的に計数する技術が求められる。

また、例えば腸管上皮のオルガノイドは、内部に空洞を有する細胞集団を構成することが知られている。これが細胞膜輸送系の評価モデルとして用いられる場合、１層の細胞層で覆われた空洞が形成されることが理想的であるが、実際には空洞の周囲に多層の細胞層が形成されることが多い。このため、腸管上皮オルガノイドの表面の細胞層の構造、具体的には細胞がどのように空洞を覆う細胞層を形成しているかを定量的に評価する技術が求められる。

このような要求に応えることが期待される技術として、本願出願人は先に特許文献１、２を開示した。特許文献１には、光干渉断層撮像（光コヒーレンストモグラフィ、Optical Coherence Tomography；ＯＣＴ）等の非侵襲の断層撮像技術により撮像した胚（受精卵）の三次元像から、栄養外胚葉と内細胞塊とを識別し分割する方法が記載されている。また、特許文献２には、ＯＣＴ撮像された胚をその重心位置を原点とする極座標により表し、各動径方向における反射光強度を輝度値に置き換えた二次元マップで視覚化することで、観察者による胚の評価作業を支援する技術が記載されている。

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

上記した従来技術は、観察対象物の三次元構造を観察者にわかりやすく提示するものであるが、観察や評価作業を支援するのに有効な定量的な情報を自動的に抽出するには至っていない。この点において、上記従来技術には改良の余地が残されている。より具体的には、胚盤胞期の胚や腸管上皮オルガノイドなど、内部に空洞を有する細胞塊の表面構造に関する定量的な情報を取得することのできる技術が求められる。

この発明は上記課題に鑑みなされたものであり、内部に空洞を有する細胞塊の画像を処理対象物とする画像処理において、ＯＣＴ撮像により得られる画像データから、細胞塊の表面構造の解析に有用な定量的情報を求めることのできる技術を提供することを目的とする。

この発明の一の態様は、内部に空洞を有する細胞塊の画像を処理対象とする画像処理方法であって、上記目的を達成するため、前記細胞塊を光干渉断層撮像して、前記細胞塊の三次元像を表す三次元画像データを取得する工程と、前記三次元画像データに基づき、前記細胞塊のうち外部空間に面する外側表面と前記空洞に面する内側表面との間の前記細胞塊の厚さを位置ごとに表す厚さ分布を取得する工程と、前記厚さ分布の取得結果に基づき、周囲の領域に比較して厚さが突出した突起部位を検出する工程とを備えている。

また、この発明の他の一の態様は、上記各工程をコンピュータに実行させるためのプログラムである。また、この発明の他の一の態様は、上記プログラムを記憶したコンピュータ読み取り可能な記録媒体である。

細胞核等の細胞の内部構造に起因して、各細胞は周縁部に比べて中央部が膨らんだ構造を有している。このため、細胞同士が互いに横方向につながって形成される細胞層では、各細胞の形状に起因する凹凸が生じる。特に胚盤胞期の胚（受精卵）のように、内部の空洞を覆うように１層の細胞により構成される細胞層では、１つ１つの細胞に対応する突起部位が当該層に現れることになる。したがって、突起部位の数が、層を構成する細胞の位置や数を表す指標となり得る。また、突起部位の間隔は個々の細胞の大きさを表す指標となり得る。

上記発明では、ＯＣＴ撮像された細胞塊の三次元画像データから、細胞塊の外側表面と、内部の空洞に臨む内側表面との間における細胞塊の厚さの分布が求められる。この厚さ分布は、空洞を覆う細胞層の厚さを位置ごとに表すことができる。細胞層の厚さは必ずしも均一ではなく、例えば１層の細胞により構成される細胞層では、個々の細胞の中央部分で厚く、周縁部ではより薄い。したがって、求められた厚さ分布から周囲に比べて厚さが大きい突起部位を検出すれば、当該部位は１つの細胞の位置を指し示すものである蓋然性が高いと言える。

このように、内部の空洞を覆う細胞層の厚さ分布を求め、そこに含まれる突起部位を検出することで、細胞層を構成する細胞の位置、数、大きさ等を定量的に表す情報を取得することが可能となる。

上記のように、本発明によれば、内部に空洞を有する細胞塊をＯＣＴ撮像して得られた画像データから、細胞塊の表面構造の解析に有用な、個々の細胞の位置や数等に関する定量的情報を求めることが可能である。

本発明の画像処理方法を実行する画像処理装置の構成例を示す図である。 本実施形態において試料となる胚の構造を模式的に示す図である。 本実施形態における画像処理を示すフローチャートである。 焦点深さと抽出される透明帯の面積との関係を示す図である。 ＸＹＺ直交座標系と極座標系との対応関係を示す図である。 栄養外胚葉の厚さの求め方を例示する図である。 栄養外胚葉の厚さを求める処理を示すフローチャートである。 二次元マッピング手法の一例を示す図である。 二次元マッピング手法の他の一例を示す図である。 三次元マッピング手法の例を示す図である。 ピーク検出処理の原理を示す図である。 ピーク検出処理を示すフローチャートである。 二次元マップ上で検出されるピーク位置の例を示す図である。 三次元マップ上で検出されるピーク位置の例を示す図である。

図１は本発明に係る画像処理方法の実行主体として好適な画像処理装置の構成例を示す図である。この画像処理装置１は、液体中に担持された試料、例えば培養液中で培養された胚（受精卵）を断層撮像し、得られた断層画像を画像処理して、試料の一の断面の構造を示す断面画像を作成する。また、複数の断面画像から試料の立体像を作成する。以下の各図における方向を統一的に示すために、図１に示すようにＸＹＺ直交座標軸を設定する。ここでＸＹ平面が水平面を表す。また、Ｚ軸が鉛直軸を表し、より詳しくは（－Ｚ）方向が鉛直下向き方向を表している。

画像処理装置１は保持部１０を備えている。保持部１０は、撮像対象物となる試料Ｓを収容する試料容器１１を水平姿勢に保持する。試料容器１１は例えばディッシュと称される平底で浅皿形状の容器である。あるいは例えば、試料容器１１は板状部材の上面に液体を担持可能な窪部（ウェル）Ｗが多数形成されたウェルプレートである。試料容器１１には培養液などの培地Ｍが注入されており、その内部に試料Ｓが担持される。

保持部１０に保持された試料容器１１の下方に、撮像部２０が配置される。撮像部２０には、被撮像物の断層画像を非接触、非破壊（非侵襲）で撮像することが可能な光干渉断層撮像（Optical Coherence Tomography；ＯＣＴ）装置が用いられる。詳しくは後述するが、ＯＣＴ装置である撮像部２０は、被撮像物への照明光を発生する光源２１と、光ファイバカプラ２２と、物体光学系２３と、参照光学系２４と、分光器２５と、光検出器２６とを備えている。

撮像部２０はさらに、光学顕微鏡撮像を行うための顕微鏡撮像ユニット２８を備えている。より具体的には、顕微鏡撮像ユニット２８は、撮像光学系２８１と、撮像素子２８２とを備えている。撮像光学系２８１は、焦点が試料容器１１内の試料Ｓに合わせられる対物レンズを含む。また撮像素子２８２としては、例えばＣＣＤ撮像素子、ＣＭＯＳセンサー等を用いることができる。顕微鏡撮像ユニット２８としては、明視野撮像または位相差撮像が可能なものが好適である。物体光学系２３と顕微鏡撮像ユニット２８とは、水平方向に移動可能な支持部材（図示省略）により支持されており、水平方向における位置を変更可能となっている。

画像処理装置１はさらに、装置の動作を制御する制御ユニット３０と、撮像ユニット２０の可動部を駆動する駆動部４０とを備えている。制御ユニット３０は、ＣＰＵ（Central Processing Unit）３１、Ａ／Ｄコンバータ３２、信号処理部３３、撮像制御部３４、インターフェース（ＩＦ）部３５、画像メモリ３６およびメモリ３７を備えている。

ＣＰＵ３１は、所定の制御プログラムを実行することで装置全体の動作を司り、ＣＰＵ３１が実行する制御プログラムや処理中に生成したデータはメモリ３７に保存される。Ａ／Ｄコンバータ３２は、撮像ユニット２０の光検出器２６および撮像素子２８２から受光光量に応じて出力される信号をデジタルデータに変換する。信号処理部３３は、Ａ／Ｄコンバータ３２から出力されるデジタルデータに基づき後述する信号処理を行って、被撮像物の画像を作成する。こうして作成された画像データは、画像メモリ３６により適宜記憶保存される。

撮像制御部３４は、撮像部２０の各部を制御して撮像処理を実行させる。具体的には、撮像制御部３４は、ＯＣＴ撮像用の物体光学系２３と光学顕微鏡撮像用の顕微鏡撮像ユニット２８とを選択的に、撮像試料Ｓを撮像視野に収める所定の撮像位置に位置決めする。物体光学系２３が撮像位置に位置決めされると、撮像制御部３４は、撮像部２０に後述するＯＣＴ撮像処理を実行させて、試料Ｓの立体構造を表す三次元画像データを取得させる。一方、顕微鏡撮像ユニット２８を撮像位置に位置決めしたときには、撮像制御部３４は、撮像光学系２８１により撮像素子２８２の受光面に結像される試料Ｓの平面像を表す二次元画像データを、光学撮像ユニット２８に取得させる。

インターフェース部３５は画像処理装置１と外部との通信を担う。具体的には、インターフェース部３５は、外部機器と通信を行うための通信機能と、ユーザからの操作入力を受け付け、また各種の情報をユーザに報知するためのユーザインターフェース機能とを有する。この目的のために、インターフェース部３５には、装置の機能選択や動作条件設定などに関する操作入力を受け付け可能な例えばキーボード、マウス、タッチパネルなどの入力デバイス３５１と、信号処理部３３により作成された断層画像や３Ｄ復元部３４により作成された立体像など各種の処理結果を表示する例えば液晶ディスプレイからなる表示部３５２とが接続されている。

撮像部２０では、例えば発光ダイオードまたはスーパールミネッセントダイオード（ＳＬＤ）などの発光素子を有する光源２１から、広帯域の波長成分を含む低コヒーレンス光ビームが出射される。細胞等の試料を撮像する目的においては、入射光を試料の内部まで到達させるために、例えば近赤外線が用いられることが好ましい。

光源２１は光ファイバカプラ２２を構成する光ファイバの１つである光ファイバ２２１に接続されており、光源２１から出射される低コヒーレンス光は、光ファイバカプラ２２により２つの光ファイバ２２２，２２４への光に分岐される。光ファイバ２２２は物体系光路を構成する。より具体的には、光ファイバ２２２の端部から出射される光は物体光学系２３に入射する。

物体光学系２３は、コリメータレンズ２３１と対物レンズ２３２とを備えている。光ファイバ２２２の端部から出射される光はコリメータレンズ２３１を介して対物レンズ２３２に入射する。対物レンズ２３２は、光源２１からの光（観察光）を試料に収束させる機能と、試料から出射される反射光を集光して光ファイバカプラ２２に向かわせる機能とを有する。図では単一の対物レンズ２３２が記載されているが、複数の光学素子が組み合わされていてもよい。被撮像物からの反射光は対物レンズ２３２、コリメータレンズ２３１を介し信号光として光ファイバ２２２に入射する。対物レンズ２３２の光軸は試料容器１１の底面に直交しており、この例では光軸方向は鉛直軸方向と一致している。

ＣＰＵ３１が撮像制御部３４に制御指令を与え、これに応じて撮像制御部３４は撮像部２０に所定方向への移動を行わせる。より具体的には、撮像制御部３４は、撮像部２０を水平方向（ＸＹ方向）および鉛直方向（Ｚ方向）に移動させる。撮像部２０の水平方向の移動により、撮像範囲が水平方向に変化する。また、撮像部２０の鉛直方向の移動により、対物レンズ２３２の光軸方向における焦点位置が、被撮像物である試料Ｓに対して変化する。

光源２１から光ファイバカプラ２２に入射した光の一部は光ファイバ２２４を介して参照光学系２４に入射する。参照光学系２４は、コリメータレンズ２４１および参照ミラー２４３を備えており、これらが光ファイバ２２４とともに参照系光路を構成する。具体的には、光ファイバ２２４の端部から出射される光がコリメータレンズ２４１を介して参照ミラー２４３に入射する。参照ミラー２４３により反射された光は参照光として光ファイバ２２４に入射する。

参照ミラー２４３は、撮像制御部３４からの制御指令によって作動する進退部材（図示省略）により支持されて、Ｙ方向に進退移動可能である。参照ミラー２４３がＹ方向、つまりコリメータレンズ２４１に対し接近および離間方向に移動することで、参照ミラー２４３により反射される参照光の光路長が調整される。

試料の表面もしくは内部の反射面で反射された反射光（信号光）と、参照ミラー２４３で反射された参照光とは光ファイバカプラ２２で混合され光ファイバ２２６を介して光検出器２６に入射する。このとき、信号光と参照光との間で位相差に起因する干渉が生じるが、干渉光の分光スペクトルは反射面の深さにより異なる。つまり、干渉光の分光スペクトルは被撮像物の深さ方向の情報を有している。したがって、干渉光を波長ごとに分光して光量を検出し、検出された干渉信号をフーリエ変換することにより、被撮像物の深さ方向における反射光強度分布を求めることができる。このような原理に基づくＯＣＴ撮像技術は、フーリエドメイン（Fourier Domain）ＯＣＴ（ＦＤ－ＯＣＴ）と称される。

この実施形態の撮像部２０は、光ファイバ２２６から光検出器２６に至る干渉光の光路上に分光器２５が設けられている。分光器２５としては、例えばプリズムを利用したもの、回折格子を利用したもの等を用いることができる。干渉光は分光器２５により波長成分ごとに分光されて光検出器２６に受光される。

光検出器２６が検出した干渉光に応じて光検出器２６から出力される干渉信号をフーリエ変換することで、試料のうち、照明光の入射位置における深さ方向、つまりＺ方向の反射光強度分布が求められる。試料容器１１に入射する光ビームをＸ方向に走査することで、ＸＺ平面と平行な平面における反射光強度分布が求められ、その結果から当該平面を断面とする試料の断層画像を作成することができる。その原理は周知であるため、詳細な説明は省略する。

また、Ｙ方向におけるビーム入射位置を多段階に変更しながら、その都度断層画像の撮像を行うことで、試料をＸＺ平面と平行な断面で断層撮像した多数の断層画像を得ることができる。Ｙ方向の走査ピッチを小さくすれば、試料の立体構造を把握するのに十分な分解能の画像データを得ることができる。これらの断層画像データから、試料の立体像に対応する三次元画像データ（いわゆるボクセルデータ）を作成することができる。

このように、この画像処理装置１は、試料容器１１に培地Ｍとともに担持された試料Ｓの画像を取得する機能を有する。画像としては、光学顕微鏡撮像により得られる二次元画像データと、ＯＣＴ撮像により得られる断層画像データおよびそれに基づく三次元画像データとを取得可能である。

以下、上記のように構成された画像処理装置１を用いて実行可能な、本発明に係る画像処理方法の一実施形態について説明する。この実施形態における画像処理は、試料Ｓとして胚盤胞期の受精卵（以下、単に「胚」という）を撮像し、その画像から、当該胚を構成する主たる構造体である透明帯、栄養外胚葉および内細胞塊に対応する領域を個別に抽出する、という内容を有するものである。

これにより得られたデータに基づき、ユーザ（具体的には医師または胚培養士）による胚の評価作業を効果的に支援することが可能となる。例えば不妊治療を目的とする胚の培養において、培養が良好に進行しているか否かを判断するための知見を得ることを目的として、本実施形態の画像処理方法を適用することができる。

図２は本実施形態において試料となる胚の構造を模式的に示す図である。既に知られているように、卵子が受精すると卵割が開始され、桑実胚と呼ばれる状態を経て胚盤胞が形成される。図２（ａ）は胚盤胞期における胚の内部構造を模式的に示すものである。胚盤胞期においては、胚Ｅは内部に胞胚腔Ｂとよばれる空洞を有している。より具体的には、卵割の進んだ細胞が胚の表面に薄い層として並び栄養外胚葉Ｔを形成し、栄養外胚葉Ｔで囲まれる内部空間が胞胚腔Ｂを形成する。また、内部空間の一部に、多数の細胞が密集した内細胞塊Ｉが形成される。

さらに、栄養外胚葉Ｔの外側表面を覆うように透明帯ＺＰが形成される。透明帯ＺＰは糖タンパク質を主体とする略均一な厚みを有する膜である。一方、多数の細胞が集まって形成される栄養外胚葉Ｔは、位置により様々な厚みを有し、透明帯ＺＰの内側表面全体に貼り付くように分布している。図２（ｂ）に示すように、栄養外胚葉Ｔは、多数の細胞Ｃが透明帯ＺＰの内面に沿って並ぶことにより形成される。個々の細胞Ｃの形状や大きさが異なることにより、栄養外胚葉Ｔは、位置により厚さが異なった層となる。なお、図２においては、栄養外胚葉Ｔの厚さのばらつきが実際よりも強調されている。

胚Ｅを構成するこれらの構造体、すなわち透明帯ＺＰ、栄養外胚葉Ｔおよび内細胞塊Ｉは、胚の評価における着目領域として重要な意味を有するものである。このため、撮像された画像からこれらの構造体に対応する領域を自動的に抽出する技術は、ユーザによる胚の評価作業を支援する上で大きな意義を有するものとなる。しかしながら、ＯＣＴ撮像により得られる断層画像あるいは三次元画像では、これらの構造体の間で輝度情報にほとんど差異がない。このため、単純な輝度の違いによる分離では、これらの構造体を精度よく分離することが困難である。

このように、ＯＣＴ画像から透明帯ＺＰ、栄養外胚葉Ｔおよび内細胞塊Ｉを自動的に分離する技術は、これまで確立されていなかった。本願出願人が先に開示した特許文献２は、このような問題に鑑み、ＯＣＴ画像から栄養外胚葉と内細胞塊とに対応する領域を分離することができるようにしたものである。ただし、この技術においても、透明帯と栄養外胚葉とが明確に分離されているとは言えなかった。具体的には、この技術により栄養外胚葉として抽出される領域が、透明帯に対応する領域まで包含している可能性がある。

本実施形態における画像処理の第１の目的は、光学顕微鏡撮像により得られた二次元画像から得られる情報を利用してＯＣＴ画像を解析することにより透明帯ＺＰと栄養外胚葉Ｔとを分離し、最終的には三次元像において透明帯ＺＰ、栄養外胚葉Ｔおよび内細胞塊Ｉが占める領域をそれぞれ個別に特定することである。また、第２の目的は、こうして分離された栄養外胚葉Ｔを構成する細胞の個数を形成することである。

光学顕微鏡撮像により得られる明視野画像または位相差画像においては、その輝度およびテクスチャの違いから透明帯ＺＰと栄養外胚葉Ｔとを区別することは比較的容易である。そして、良好に培養された胚盤胞期の胚では、透明帯ＺＰは概ね均一な厚みを有している。言い換えれば、胚のうち、その表面から一定の深さまでの領域は透明帯ＺＰが占めていると考えることができる。これらのことから、光学顕微鏡撮像により得られた二次元画像データから透明帯ＺＰの厚さを特定し、ＯＣＴ撮像により得られた三次元画像データから、胚Ｅの表面から見て透明帯ＺＰの厚さに相当する深さまでの領域を特定することで、透明帯ＺＰに対応する領域を栄養外胚葉Ｔとは分離して抽出することが可能である。また、内細胞塊Ｉの分離については、例えば特許文献２に記載された方法を適用可能である。これらにより、第１の目的を達成することができる。

こうして他の構造体から分離された栄養外胚葉Ｔは、図２（ｂ）に示したように、個々の細胞Ｃの形状に起因する凹凸を有している。より具体的には、栄養外胚葉Ｔの胞胚腔Ｂに臨む表面には、個々の細胞Ｃの形状に起因する突起が生じている。したがって、このような突起部位を検出し、その数を計数すれば、栄養外胚葉Ｔを構成する細胞Ｃの個数を知ることができる。このようにして第２の目的が達成される。

図３は本実施形態における画像処理を示すフローチャートである。この処理は、ＣＰＵ３１が予め用意された制御プログラムを実行し、装置各部に所定の動作を行わせることにより実現される。評価対象となる胚が収容された試料容器１１が保持部１０にセットされると（ステップＳ１０１）、当該胚を撮像対象物として、撮像部２０による光学顕微鏡撮像およびＯＣＴ撮像が実行される。

図３において、ステップＳ１０２ないしＳ１０５は、光学顕微鏡撮像およびこれにより得られる二次元画像データを用いたデータ処理を表す。一方、ステップＳ１１１およびＳ１１２は、ＯＣＴ撮像およびこれにより得られる三次元画像データを用いた画像処理を表す。これら２つの処理は、時系列的に前後して行われてもよく、この場合いずれの処理を先に行うかは任意である。また、これらの処理を併行して行うことで処理時間の短縮が図られてもよい。

まずステップＳ１０２ないしＳ１０５の処理の概要について説明する。この処理では、胚を光学顕微鏡撮像することで胚の二次元画像データを取得し、画像データに基づき透明帯の厚さを算出する。具体的には、顕微鏡撮像ユニット２８が撮像位置に位置決めされ、焦点位置が深さ方向（Ｚ方向）に多段階に変更設定されるとともにその都度撮像が行われる。これにより、焦点深さを互いに異ならせた複数の画像セット、いわゆるＺスタック画像が取得される（ステップＳ１０２）。

これらの画像の各々から、透明帯ＺＰに対応する領域が抽出される（ステップＳ１０３）。そして、複数の画像のうち、透明帯ＺＰに最もよく合焦した１枚が選出され（ステップＳ１０４）、選出された画像に基づき、透明帯ＺＰの厚さが算出される（ステップＳ１０５）。これらの処理、すなわち透明帯ＺＰに対応する領域を抽出する処理（ステップＳ１０３）、最も合焦した画像を選出する処理（ステップＳ１０４）および透明帯ＺＰの厚さを算出する処理（ステップＳ１０５）の詳細については後述する。ここまでの処理により、評価対象の胚Ｅにおける透明帯ＺＰの厚さが既知となっている。

一方、ステップＳ１１１では、物体光学系２３が撮像位置に位置決めされ、撮像位置を走査させながら胚Ｅが断層撮像されることにより、胚Ｅの三次元画像データが取得される。三次元画像データで示される各部の輝度が所定の閾値で二値化されることにより、胚Ｅの三次元像は、比較的高密度な何らかの構造体で占められて高輝度となる領域（以下、「構造領域」と称する）と、これより低密度で低輝度の領域とに分離される。例えば培地ＭはＯＣＴ画像において低輝度である。

構造領域を占める構造体としては、前述の通り、透明帯ＺＰ、栄養外胚葉Ｔおよび内細胞塊Ｉを含み得る。ステップＳ１１３、Ｓ１１４ではこれらが相互に分離される。具体的には、ステップＳ１１３において、透明帯ＺＰに対応する領域が他の構造領域とは分離される。そして、ステップＳ１１４において、栄養外胚葉Ｔと内細胞塊Ｉとが分離される。

ここまでの処理（ステップＳ１０１～Ｓ１１４）により、胚Ｅの三次元像から透明帯ＺＰ、栄養外胚葉Ｔおよび内細胞塊Ｉが占める領域をそれぞれ特定する、という第１の目的が達成される。続くステップＳ１１５～Ｓ１１８の処理は、分離された栄養外胚葉Ｔを構成する細胞の個数を計数するという第２の目的のためのものである。この処理については後で詳しく説明する。

以下では、上記処理の各工程（ステップＳ１０３～Ｓ１０５、Ｓ１１１）を実行するための各要素技術について順に分説する。なお、ステップＳ１０２のＺスタック画像を取得するための処理は周知であるため、説明を省略する。また、ステップＳ１１４の処理についても、特許文献２に記載された技術を適用可能であるため説明を省略する。

ステップＳ１０３では、胚Ｅを光学顕微鏡撮像して得られた二次元画像データから、透明帯ＺＰに対応する領域が抽出される。この処理は適宜の画像処理技術を用いて実行することが可能である。例えば、画像内から特定の特徴を有する領域を抽出するパターン認識技術を適用することができる。具体的には、予め取得された透明帯の画像を教師画像とする教師付き学習により分類モデルを構築し、これを用いて評価対象の胚Ｅの光学顕微鏡画像を領域分割することにより、画像から透明帯ＺＰに対応する領域を抽出することができる。

この実施形態では、領域分割処理の一例として公知のセマンティックセグメンテーション法を用いる。セマンティックセグメンテーション法は、ディープラーニングアルゴリズムにより予め構築された分類モデルを用いて、画像内の各画素にラベル付けを行う技術である。本実施形態では、次のようにしてこの方法を利用することができる。

まず、透明帯が良好な画像品質で撮像された胚の光学顕微鏡画像を用意し、当該画像中の透明帯に対応する領域の各画素に、その旨を表すラベル付けを行う。そして、元の顕微鏡画像を入力データ、ラベル画像を正解データとしてディープラーニングを実行することで、分類モデルが構築される。こうして予め構築された分類モデルに未知の画像を入力データとして与えることで、当該画像のうち透明帯に対応する領域にその旨のラベルが付された出力画像を得ることができる。このような領域を出力画像から抽出することで、結果的に透明帯を抽出することが可能である。

分類モデルを構築するための具体的手法の一例を説明する。この処理は、画像を表示する機能およびユーザからの操作入力を受け付ける機能を有する各種のコンピュータ装置により実行可能である。例えば画像処理装置１、あるいはパーソナルコンピュータ等の汎用コンピュータ装置によって実行することが可能である。

最初に、透明帯に合焦した（ピントが合った）状態で予め撮像された、胚の光学顕微鏡画像を表示する。画像処理装置１では、表示部３５２に当該画像を表示させることができる。こうして表示された画像に対し、ピントの合った透明帯に相当する領域を指定するためのユーザからの教示入力を受け付ける。この場合のユーザは、胚の画像に対する十分な知見を有する熟練者であることが望ましい。また、画像処理装置１が用いられる場合、入力デバイス３５１を介して教示入力を受け付けることができる。

透明帯として指定された領域にはその旨を示すラベルが付与される。こうしてラベル付けされた画像を正解データとし、元の画像を入力データとしてディープラーニングを実行することで、画像から透明帯を抽出するための分類モデルが構築される。必要に応じ、透明帯以外のラベルが用いられてもよい。

この分類モデルは、透明帯にピントの合った画像を入力画像として構築されたものである。そのため、これを未知のテスト画像に適用してセマンティックセグメンテーション法を実行すると、テスト画像から透明帯としての特徴を強く有する領域が抽出されることになる。こうして透明帯が抽出されると、その厚さを求めることが可能となる。厚さを精度よく求めるためには、画像内のできるだけ広い領域で、透明帯に合焦していることが望ましい。すなわち、透明帯の厚さの算出は、セマンティックセグメンテーション法で抽出される領域の面積ができるだけ広い画像に基づき行われることが望ましいと言える。なお面積については画素数によって表すことが可能である。

一方、評価対象の胚Ｅは立体的構造を有しているため、焦点深さを適宜に定めて撮像した画像では、透明帯ＺＰへの合焦状態が必ずしも良好でない場合があり得る。そこで、深さ方向（Ｚ方向）に焦点位置を異ならせて取得されたＺスタック画像から、透明帯ＺＰに相当するものとして抽出される領域の面積が最も大きい１つの画像を選出し、この画像を透明帯ＺＰの厚さ算出に用いる。

ステップＳ１０４では、Ｚスタック画像から透明帯ＺＰに最も合焦した画像が選出される。この実施形態におけるセマンティックセグメンテーション法では、画像からピントが合った透明帯の領域が抽出されるから、この領域の面積が最大である画像が、透明帯ＺＰに最も合焦した画像である蓋然性が高いと言える。そのような条件に該当する画像を選出すればよい。

図４は焦点深さと抽出される透明帯の面積との関係を示す図である。図４（ａ）の上部には、深さ方向における種々の焦点位置で撮像される光学顕微鏡画像（Ｚスタック画像）と、それらの画像からセマンティックセグメンテーション法により抽出される領域（透明帯ＺＰに相当する領域）との関係が模式的に示されている。また、図４（ａ）の下部のグラフは、撮像時の焦点位置と抽出された領域の面積（画素数）との関係が例示されている。これらに示されるように、Ｚスタック画像のうち最もよく合焦し透明帯ＺＰが鮮明である画像Ｉａにおいて、抽出される領域の面積が最大になる。言い換えれば、Ｚスタック画像のうち抽出された面積が最大である画像に対応する焦点位置を合焦位置とみなすことができる。したがって、この画像Ｉａを、透明帯ＺＰの厚さ算出の基となる画像として選出すればよい。

ただし、例えば撮像時の振動等に起因して、画像にブレが生じ、結果として透明帯の見かけ上の面積が実際より大きく画像に現れる場合があり得る。図４（ｂ）はこのような状況を表したものである。例えば画像Ｉｂにおいて、振動に起因する画像のブレが生じているとする。そうすると、透明帯に相当するとして抽出される面積は見かけ上大きくなり、この位置が合焦位置と誤判定されるおそれがある。

この実施形態では、抽出された領域の周縁部を挟む画素の間の輝度差を用いることで、この問題の解消が図られている。すなわち、よく合焦した画像では、透明帯に対応する領域とその周囲領域との境界が明りょうであり、したがってそれらの領域の輝度の間にも明確なコントラストがあると考えられる。一方、合焦していない画像では、これらの領域間の境界が明確でなく、したがってコントラストも鮮明でない。

このことから、単に抽出された領域の面積（画素数）の大小だけで評価するのに代えて、当該領域のエッジ部分における輝度変化の大きさ、すなわちシャープネスも反映させた評価値を導入することで、上記のような誤判定を低減することができると考えられる。このようなエッジ変化量を定量化する方法としては種々のものがあり、それらを適宜選択して適用可能である。

この実施形態ではその一例として、抽出された領域の面積に、エッジ部分における輝度変化の大きさを反映した係数を乗じた値を評価値とする。この係数としては、例えばエッジを挟んだ画素間における輝度の差を二乗した値を用いることができる。より具体的には、抽出された領域のうちエッジに隣接する全ての画素の輝度の平均値と、当該エッジの外側でエッジに隣接する全ての画素の起動の平均値との差を求め、それを二乗することで上記係数が求められる。

このようにすると、図４（ｂ）下部のグラフに示されるように、画像Ｉｂにおけるズレの影響による抽出面積の増大が合焦位置の誤判定の原因となるリスクを低減させることができる。なお、ここでは係数を正の値とするために輝度差を二乗しているが、これに代えて輝度差の絶対値を係数として用いてもよい。

透明帯に合焦した画像を選出するための処理は、以下のように構成することができる。この処理は、図３のステップＳ１０４に相当するものである。Ｚスタック画像を構成する各画像に対しては、ステップＳ１０３において透明帯に相当する領域が抽出されている。こうして抽出された領域の面積を求めるため、当該領域に属する画素の数を画像ごとに計数する。そして、抽出された領域のエッジの内外においてそれぞれ隣接する画素の平均輝度を求め、その差を算出する。これらの値に基づき、各画像の合焦度合いを示す評価値を算出する。具体的には、抽出領域の画素数に、エッジ内外での輝度差の二乗で表される係数を乗じることで、評価値を算出する。こうして求められる評価値が最大の画像が、透明帯に最も合焦した画像として選出される。

次に、Ｚスタック画像から選出された１つの光学顕微鏡画像から透明帯ＺＰの厚さを算出する、ステップＳ１０５の処理内容について説明する。選出された画像については、ステップＳ１０３において透明帯ＺＰに相当する領域の抽出が行われる。良好に培養された胚Ｅにおいては、概ね一定の幅を有する円環状の領域が、透明帯ＺＰに相当する領域として抽出されていると考えられる。透明帯ＺＰに合焦した画像では、この円環の幅、つまり円環の内側のエッジと外側のエッジとの距離が、透明帯ＺＰの厚さに対応している。

円環の幅を求める方法としては種々のものが考えられるが、簡便なものとしては、ＯｐｅｎＣＶ（Open Source Computer Vision）ライブラリに装備されているDistance Transform関数を利用した方法がある。円環の内側エッジ上の１つの画素を注目画素としてDistance Transform関数を適用することで、当該画素から最も近い、円環の外側エッジ上の画素までの距離を特定することができる。この距離が当該位置における円環の幅、つまり透明帯ＺＰの厚さを表す。これとは逆に、円環の外側エッジ上の画素を注目画素とし、これから内側エッジまでの最短距離を求めても等価である。こうして円環上の各位置で求められる幅の平均値を、透明帯ＺＰの厚さの平均値とすることができる。以下ではこの厚さの平均値を符号Ｔａにより表すこととする。

以上、ステップＳ１０３～Ｓ１０５の詳細な処理内容について説明した。次に、ステップＳ１１３における、透明帯の分離処理について説明する。ここでは、ＯＣＴ撮像により得られた胚Ｅの三次元像のうち、表面から深さＴａまでの範囲を透明帯ＺＰとみなす。したがって、三次元像からこの範囲のみを抽出することで、透明帯ＺＰに対応する構造体のみを取り出すことができる。一方、この範囲の構造体を消去することで、透明帯ＺＰ以外の構造体、すなわち栄養外胚葉Ｔと内細胞塊Ｉとを取り出すことができる。

透明帯ＺＰに相当する領域が消去された三次元像に対しては、特許文献２に記載されているように、残された構造体に対し例えばLocal Thickness演算を利用した領域分割処理を実行することで、栄養外胚葉Ｔに相当する領域と、内細胞塊Ｉに相当する領域とを分離することができる。このようにして、胚Ｅの三次元像から、透明帯ＺＰ、栄養外胚葉Ｔおよび内細胞塊Ｉがそれぞれ占める領域を個別に抽出することができる。

次に、本実施形態の第２の目的を達成する、すなわち栄養外胚葉Ｔを構成する細胞数を計数するための処理（ステップＳ１１５～Ｓ１１９）について説明する。この処理の概略の流れは、栄養外胚葉Ｔの各位置における厚さ分布を算出し、求められた厚さ分布におけるピーク位置を検出し、有意なピークが個々の細胞に対応するものとしてその個数を計数する、というものである。以下、各工程について順に説明する。

ステップＳ１１５では、ＸＹＺ直交座標系で表されている三次元画像データが極座標表現に変換される。胚Ｅは概ね球形で内部が空洞であるため、栄養外胚葉Ｔは球殻に近い形状を有している。このような構造をより端的に表現するためには、胚Ｅの中心を原点とする極座標（球座標）で各位置を表す方が好ましい。そこで、図５に示すように、ＸＹＺ直交座標系から、動径ｒおよび２つの偏角θ、φを座標変数とするｒθφ極座標系への座標変換を行う。

図５はＸＹＺ直交座標系と極座標系との対応関係を示す図である。よく知られているように、直交座標系における点Ｐの座標（ｘ，ｙ，ｚ）と球座標系における点Ｐの座標（ｒ，θ，φ）との間には、原点Ｏが共通であれば下式：
ｘ＝ｒ・ｓｉｎθ・ｃｏｓφ
ｙ＝ｒ・ｓｉｎθ・ｓｉｎφ
ｚ＝ｒ・ｃｏｓθ
の関係がある。

具体的には、三次元画像データから胚Ｅの中心位置を特定し、この中心位置を球座標系の原点Ｏとする。この原点Ｏは、原信号データにおけるＸＹＺ直交座標系の原点とは一致する必要はない。そして、適宜の変換処理により直交座標から球座標への座標変換を行う。このように座標変換を行うことにより、ＸＹＺ直交座標系で特定される三次元空間内の各位置を、極座標表現することができる。

「胚の中心」については、三次元画像データに基づき例えば次のようにして求めることができる。三次元画像データで表される胚Ｅの三次元像の表面が球面とみなせる場合には、画像における当該球の重心を胚の中心とすることができる。三次元画像処理において、中実のオブジェクトの重心を求める方法は公知であり、そのような方法を適用することが可能である。また、例えば胚の表面を近似的に表す球面または回転楕円面を特定し、該近似面の中心を胚の中心としてもよい。

ステップＳ１１６では、こうして極座標変換された三次元画像データから、栄養外胚葉Ｔの厚さ分布が求められる。具体的には、極座標空間内の１つの動径方向における栄養外胚葉Ｔの厚さを算出し、これを種々の動径方向について行うことにより、栄養外胚葉Ｔの厚さ分布を求めることができる。

図６は栄養外胚葉の厚さの求め方を例示する図である。図６（ａ）に示す例では、原点Ｏから延びる一の動径ｒ１を考え、当該動径ｒ１方向における栄養外胚葉Ｔの内側表面つまり胞胚腔Ｂに臨む面Ｓｉと、外側表面つまり透明帯ＺＰに臨む面Ｓｏとの距離Ｔ１を、当該動径方向における栄養外胚葉Ｔの厚さとする。

また、図６（ｂ）に示す例では、動径ｒ１が栄養外胚葉Ｔの内側表面Ｓｉに交わる点Ｐ１に着目し、当該点Ｐ１と栄養外胚葉Ｔの外側表面Ｓｏとの最短距離を、当該動径ｒ１方向における栄養外胚葉Ｔの厚さとする。あるいは、動径ｒ１が栄養外胚葉Ｔの外側表面Ｓｏに交わる点Ｐ２に着目し、当該点Ｐ２と内側表面Ｓｉとの最短距離を、当該動径ｒ１方向における栄養外胚葉Ｔの厚さとしてもよい。最短距離の算出には、前記した透明帯ＺＰの厚さ算出処理と同様に、ＯｐｅｎＣＶライブラリ中のDistance Transform関数を利用することが可能である。

図７は栄養外胚葉の厚さを求める処理を示すフローチャートである。最初に、一の動径方向が選択される（ステップＳ２０１）。具体的には、偏角θ，φの値が適宜に仮設定されることにより、１つの動径方向が特定される。次に、極座標空間内で当該動径方向における栄養外胚葉Ｔの内側表面Ｓｉの位置（図６（ｂ）に示す点Ｐ１の位置）、または外側表面Ｓｏの位置（図６（ｂ）に示す点Ｐ２の位置）に対応する画素が、注目画素として特定される（ステップＳ２０２）。後の処理のために、このときの原点Ｏから注目画素までの距離を、偏角θ，φの組で特定される動径方向と関連付けて記憶させておくことが望ましい。この情報を関数Ｒ（θ，φ）により表すこととする。

続いて、当該動径方向における栄養外胚葉Ｔの内側表面Ｓｉと外側表面Ｓｏとの距離、すなわち栄養外胚葉Ｔの厚さが求められる（ステップＳ２０３）。厚さは図６（ａ）、図６（ｂ）のいずれかに示した方法で求めることができる。

このようにして求められる栄養外胚葉Ｔの厚さは１つの動径方向と関連付けられており、したがって偏角θ，φの関数として表すことができる。以下では、偏角θ，φの組により特定される１つの動径方向における栄養外胚葉Ｔの厚さをＴ（θ，φ）と表すこととする。求められた厚さＴ（θ，φ）はメモリ３７に記憶保存される（ステップＳ２０４）。種々の動径方向について上記処理を繰り返すことにより（ステップＳ２０５）、各方向における栄養外胚葉Ｔの厚さを表す厚さ分布を取得することができる。

ステップＳ１１７では、こうして求められた栄養外胚葉Ｔの厚さの空間分布を表すマップを作成する。厚さ分布マップとしては、二次元マップおよび三次元マップが考えられる。また二次元マップとしては、全ての方向を１つのマップに収める方法と、極座標空間をいくつかの空間に分割し、それぞれを１つずつのマップに収める方法とが考えられる。以下、これらのマッピング方法の具体的な作成例につき説明する。

図８は二次元マッピング手法の一例を示す図である。極座標で表現される各方向を二次元平面上にマッピングするには、例えば擬円筒図法を用いることが可能である。すなわち、図８（ａ）に示すように、偏角θ、φをそれぞれ緯度、経度に相当するものとして、略球体である地球表面の地形を平面地図上に表現する場合と同様の図法を用いることができる。このような作図法はいくつか知られており、目的に応じて適宜選択することができる。

ここでは、偏角θを緯度、偏角φを経度として扱い、実体における面積を図上で維持することのできる地図図法であるサンソン図法が用いられている。偏角θ，φの組み合わせで特定されるマップ上の各点には、各動径方向において算出された栄養外胚葉Ｔの厚さＴ（θ，φ）に対応する輝度値を有する画素が配置される。このようにして各方向での厚さをマップ上に表すことで、厚さ分布を可視化することができる。図８（ｂ）は、各点での栄養外胚葉Ｔの厚さを輝度で表したマップを、直接輝度により表現するのに代えて等輝度線を用いて表したものである。単に各方向での厚さを等高線で表したものとも言える。この他、栄養外胚葉Ｔの厚さを色や濃淡によって表すことも可能である。

なお、図８（ｂ）においてハッチングを付した右上の領域は、内細胞塊Ｉが存在する領域を表している。内細胞塊Ｉと透明帯ＺＰとの間にも栄養外胚葉Ｔは存在するが、内細胞塊Ｉに接している表面の形状は必ずしも個々の細胞の形状を表さない。このため、この実施形態では、内細胞塊Ｉが存在する領域については厚さ分布の導出対象から除外している。

地図作成の場合においても問題となっているように、略球体の表面を二次元マップによって表す場合、全ての情報を正確に表現することができない。例えば正積擬円筒図法の一種であるサンソン図法の場合、地表の面積を正確に表現できるが、距離や方向は必ずしも正しく表現できない。本実施形態のマッピングの場合、特にマップ平面の周縁部において歪みが大きくなる。

図９は二次元マッピング手法の他の一例を示す図である。図９（ａ）に、極座標空間で略球面として表される栄養外胚葉Ｔの三次元像を、原点Ｏを通り互いに垂直な２つの平面によって４つに分割する。そして、分割された各区画Ｔａ～Ｔｄをそれぞれ１つのマップに表す。このとき、区画Ｔａに対応するマップでは、区分された球面の中心、より具体的には、原点Ｏを通り、かつ２つの平面のそれぞれに対し４５度の角度で交わる直線と球面との交点Ｃａが、マップにおいて概ね中心となるようにする。他の区画Ｔｂ，Ｔｃ，Ｔｄについても同様である。

このようにすると、図９（ｂ）に示すように、極座標空間で４分割された栄養外胚葉Ｔが、マップ平面の中央部のみを使って表現されることになる。このため、マップは４面となるものの、各マップでは、歪みが大きくなる周縁部を使わずに栄養外胚葉Ｔの全体における厚さ分布を精度よく表現することが可能となる。

なお、現実的には、栄養外胚葉Ｔに対応する球面の分割を考える際に各区画を部分的にオーバーラップさせておき、マッピングの際には、図９（ｂ）に点線で示すように区画Ｔａ内のみならずその外側の一定範囲まで厚さ分布が示されていることがより好ましい。マッピングの便宜上の理由により分割されるものの、実際の栄養外胚葉Ｔは各区画よりも外側まで連続しており、区画外の一定範囲までマッピングを行っておくことで、各区画の外縁における厚さ分布の連続性を表現することができるからである。

この意味においては、分割を行うのに代えて、マップの中心となる位置を異ならせた複数の二次元マップを作成し、それらの中央部分のみを使用するようにすることも考えられる。また、分割数も上記の４に限定されず任意である。

上記のような二次元マッピング手法は、例えば画面上あるいは印刷紙面において全体の厚さ分布を俯瞰するような用途に適しているが、その周縁部において正確性が低下することは避けられない。このため、精度の面では三次元マップの方が優れている。三次元マップにおいては、極座標空間における全ての方向について、同程度の精度を確保することが可能である。

三次元マップを画面表示するのに際しては、例えば、ある１つの視線方向に投影されたマップを一時的に表示させておき、ユーザ操作に応じて視線方向を変化させるようにすればユーザの利便性を向上させることが可能である。例えばユーザ操作に連動してリアルタイムで視線方向が変化するアニメーション表示とすれば、ユーザはあたかも眼前の胚Ｅを直接観察しているかのような感覚で評価作業を行うことが可能となる。

栄養外胚葉Ｔについての三次元マッピングの手法としては、栄養外胚葉Ｔの概略形状を表す近似球面（または回転楕円面）に各方向における厚さ分布を表す情報を付与する方法と、実際の栄養外胚葉Ｔの三次元形状に相当する曲面に厚さ分布の情報を付与する方法とがある。さらに後者の場合、曲面の形状を、栄養外胚葉Ｔの内側表面Ｓｉに相当するものとするか、外側表面Ｓｏに相当するものとするかの選択肢がある。

栄養外胚葉Ｔの外側表面Ｓｏは透明帯ＺＰに接しており、したがってその表面形状は比較的滑らかである。この意味において、栄養外胚葉Ｔの外側表面Ｓｏに相当する曲面は、実質的には球面とあまり変わらない。一方、栄養外胚葉Ｔの内側表面Ｓｉは、栄養外胚葉Ｔを構成する個々の細胞の形状に対応する凹凸を有している。このため、栄養外胚葉Ｔの内側表面Ｓｉに相当する曲面にマッピングを行う方が、実際の凹凸形状をイメージしやすいという利点がある。

なお、マップ表面で厚さを可視化する方法としては、二次元マッピングの場合と同様に、厚さを輝度に変換して表示する方法、等高線（等輝度線）や色分け、濃淡等による表現方法等を適用可能である。以下では、球面にマッピングを行うケースと、栄養外胚葉Ｔの内側表面Ｓｉに相当する曲面にマッピングを行うケースとを例示する。マッピング方法としては等高線を用いたものとする。

図１０は三次元マッピング手法の例を示す図である。図１０（ａ）は球面にマッピングを行った例、図１０（ｂ）は栄養外胚葉Ｔの内側表面Ｓｉに相当する曲面にマッピングを行った例をそれぞれ示す。ここでは三次元マップが二次元平面に投影された状態で示されているが、実際には三次元曲面上でのマッピングがなされている。画面への表示出力の場合、視線方向を種々に異ならせた表示を行うことが可能である。なお、これらのマップにおいても、ハッチングで示される領域は内細胞塊Ｉに対応している。

例えば図１０（ａ）の三次元マップは、極座標空間中の原点Ｏを中心とする球面上の各位置に、当該位置を通る動径の方向における厚さ算出結果に対応する輝度の画素を配置することにより作成可能である。すなわち、２つの偏角（θ，φ）で表される動径と球面との交点に、厚さＴ（θ，φ）に対応する輝度値の画素を配置すればよい。

三次元マップでは二次元マップで生じる歪みの問題が解消される。二次元マップでは歪みの問題とは別に、動径の長さ、つまり原点Ｏからの距離の情報がマップに反映されないという問題がある。特に試料Ｓの形状が球面から大きく乖離したものである場合には、マップに表示される凹凸が実際の形状をよく反映したものとならない場合があり得る。三次元マップであれば、例えば試料の形状により近い回転楕円体を近似曲面とすることで、このような乖離を低減することが可能である。

特に図１０（ｂ）に示すように、実際の形状を反映させた曲面にマッピングを行うことで、より現実に即したマップを作成することができる。マッピングのために必要な情報は、各動径方向における栄養外胚葉Ｔの内側表面Ｓｉ（または外側表面Ｓｏ）の位置（動径の長さ）と、当該方向における栄養外胚葉Ｔの厚さとである。これらの情報は、栄養外胚葉Ｔの厚さを算出する工程（図３のステップＳ１１６）においてそれぞれＲ（θ，φ）、Ｔ（θ，φ）として既に求められている。したがって、極座標空間において座標（Ｒ（θ，φ），θ，φ）で表される位置に、厚さＴ（θ，φ）に対応する輝度値を有する画素を配置することで、三次元マップが作成される。

なお、図１０（ｂ）に示す三次元マップを二次元の表示画面に投影表示する場合、マップ曲面の凹凸形状と各位置での栄養外胚葉Ｔの厚さとを識別可能に表示する必要がある。特に投影された二次元図形の外周部以外の曲面の形状については、表示画像からは読み取りにくいものとなる可能性がある。この問題を回避するためには、曲面の形状と各点での厚さとを異なる方法で表現する必要がある。例えば曲面の形状を等高線によって表す一方、各部での厚さを濃淡や色分けにより表すようにすればよい。

図３に戻って、このようにしてマッピングされた厚さ分布から栄養外胚葉Ｔを構成する細胞の数を計数する方法についての説明を続ける。上記のようにして可視化される栄養外胚葉Ｔの表面の凹凸は個々の細胞の形状に起因するものである。このことから、例えば厚さが輝度で表現されている場合、マップ上に散在している周囲よりも高輝度の領域は、栄養外胚葉Ｔ内で１つ１つの細胞に対応する突起部位に対応しているものと考えられる。したがって、このような周囲より厚い、つまり高輝度の領域を抽出すれば、個々の細胞がどこにあるかを推定することができ、また栄養外胚葉Ｔを構成する細胞の個数を計数することができる。

このことから、栄養外胚葉Ｔの厚さプロファイルに現れるピークは、栄養外胚葉Ｔを構成する個々の細胞に対応する突起部位を示すものである可能性が高い。そこで、ステップＳ１１８では厚さプロファイルにおける有意なピークを見出すことで、間接的に突起部位を検出する。

図１１はピーク検出処理の原理を示す図である。また、図１２はピーク検出処理を示すフローチャートである。図１１（ａ）に実線で示すように、マップ上のある方向に沿って当該方向の位置と輝度（すなわち厚さ）との関係（厚さプロファイル）をプロットすると、位置ごとに輝度の変動があり、周囲に対して輝度が突出して高いピークが現れる。これらは細胞の位置に対応するピークの候補となるものである。

まず、このように周囲よりも輝度の高い領域を抽出する。具体的には、実線で示される厚さプロファイルに対し最大値フィルタ処理を実行する（ステップＳ３０１）。ウィンドウサイズを適宜に設定して最大値フィルタ処理を実行することで、図１１（ａ）に点線で示すように、プロファイルのピーク幅が膨張する。最大値処理によって、ピーク以外の位置では輝度値が上昇するが、ピーク位置では輝度値に変化はない。そこで、このように輝度値の変化のない位置（図に破線で示す）を、ピーク候補として抽出する（ステップＳ３０２）。

図１１（ｂ）は抽出されたピーク候補を示している。ピークにおける輝度が極端に高い、あるいは極端に低い場合、当該ピークは、画像ノイズや細胞に対応しないピークなど検出対象とは異なるものである可能性がある。そこで、輝度値に対し閾値Ｌth1、Ｌth2を設定し、ピーク候補のうち輝度値が上限値を規定する閾値Ｌth1を超えるもの、および輝度値が下限値を規定する閾値Ｌth2を下回るものについてはこれを消去する（ステップＳ３０３）。図１２ではこの処理を「上限・下限処理」と記載している。図１１（ｂ）において、黒丸印はこの処理により消去されるピーク候補を示している。また白丸印は消去されずに残るピーク候補を示している。

栄養外胚葉Ｔを構成する個々の細胞を特定するという本実施形態の目的においては、閾値Ｌth1、Ｌth2については、栄養外胚葉Ｔを構成する細胞として妥当なサイズの上限値および下限値に対応して設定することができる。これにより、１つの細胞と見なせないような極端に大きな、あるいは小さなピーク候補を除外することができ、細胞数の計数結果に及ぼす誤差を低減させることができる。なお、ここではピーク輝度値に上限値および下限値を設定しているが、いずれか一方であってもよい。また、大きな誤差要因とならないのであれば、上限・下限処理を省いてもよい。

一方、栄養外胚葉Ｔの面に沿った方向においても妥当な細胞のサイズがある。言い換えれば、ピーク候補が個々の細胞を適切に表すものであるならば、それらのピーク間の距離は所定の範囲内に収まっているはずである。そこで、残るピーク候補のうち互いに隣接するものの間の距離を算出し（ステップＳ３０４）、極端に近接したピーク候補がある場合には、それらの一方をノイズとみなして消去する（ステップＳ３０５）。より小さなピークを消去するのが現実的である。図１１（ｃ）は、隣接ピークとの距離が予め設定された最小距離Ｄminに満たないピーク候補を消去した結果の例を示している。

このようにして、図１１（ｃ）において白丸印で示されるピークが最終的に残る。これらのピーク位置が、個々の細胞に対応する突起部位の位置を表しているものと考えられる。つまり、ここまでの処理により、胞胚腔Ｂを取り囲む栄養外胚葉Ｔに存在する突起部位が検出されたことになる。

なお、ここでは概念を理解しやすくするために一次元の厚さプロファイルを用いて原理説明を行ったが、実際の厚さプロファイルは、二次元マップを用いる事例では二次元、三次元マップを用いる事例では三次元である。このため、二次元マップを用いる事例では、ステップＳ３０１における最大値処理は二次元フィルタリングとなる。各方向に一定のウィンドウサイズを適用するため、この場合には円形フィルタが用いられる。また、ピーク間の距離についても、二次元マップ平面上での距離となる。

二次元マップでは、動径の長さの情報が反映されておらず、また周縁部に近づくほど歪みが大きくなる。このため厳密には、円形のフィルタウィンドウで理想的なフィルタリング処理が実行されるとは言えない。しかしながら、個々の細胞位置に対応するピークを検出するという目的においては、実用的には十分な精度を得ることが可能である。特に、図９に示すように対象物を複数に分割してマッピングする方法では、歪みが大きくなるマップの周縁部を使用しないため、このような要因による誤差を十分に抑えることが可能である。図９（ｂ）に点線で示す範囲については、少なくともフィルタウィンドウのサイズよりも大きくしておくことで、フィルタ処理によって誤差が生じることは避けられる。

また、厚さプロファイルが三次元マップで示される場合、フィルタウィンドウは三次元の球形となる。すなわち、最大値フィルタ処理は球形フィルタリングとなる。また、ピーク間の距離も三次元空間内での距離となる。

各位置での栄養外胚葉Ｔの厚さが球面に投影された三次元マップ（図１０（ａ））では、原点Ｏから注目画素までの動径の長さの情報が処理に反映されない。このことは処理を簡単にする効果を有するが、フィルタ処理における誤差要因ともなり得る。一方、栄養外胚葉Ｔの表面形状に対応させた曲面へ投影される三次元マップ（図１０（ｂ））では、注目画素の位置には原点Ｏからの距離が反映されているため、当該注目画素を中心とする球形フィルタ処理により、誤差を排除し適切にピーク検出を行うことが可能である。

図１３は二次元マップ上で検出されるピーク位置の例を示す図である。図１３（ａ）に点線で示すように栄養外胚葉Ｔの厚さが周囲に比べて局所的に突出している領域のうち、上記した条件に該当するものが有意なピーク（黒丸印）として検出され、これらのピークは、栄養外胚葉Ｔを構成する細胞の１つ１つに対応していると考えられる。図１３（ｂ）は検出された有意なピークのみを二次元マップにプロットしたものである。このように、個々の細胞に対応する有意なピークが検出されると、それらの個数を計数することで突起部位の数が求められ（ステップＳ１１９）、間接的に栄養外胚葉Ｔを構成する細胞の数を知ることができる。

内細胞塊Ｉが占める領域においては、栄養外胚葉Ｔを構成する細胞の分布を知ることは困難であるが、他の領域と同程度の密度で細胞が分布していると仮定することで、胚Ｅ全体としての細胞の個数を推定することが可能である。すなわち、計数された細胞の個数を内細胞塊Ｉと接する部分を除く栄養外胚葉Ｔの表面積で除することにより細胞密度を求め、内細胞塊Ｉと接する部分を含む栄養外胚葉Ｔの表面積に細胞密度を乗じることで、栄養外胚葉Ｔの全体を構成する細胞の総数を見積もることが可能である。

図１４は三次元マップ上で検出されるピーク位置の例を示す図である。三次元マップにおいても上記と同様に、栄養外胚葉Ｔの厚さが周囲に比べて局所的に突出している領域のうち、条件に該当するものが細胞に対応する有意なピーク（黒丸印）として検出される。これをマップ曲面の全体において計数することで、栄養外胚葉Ｔを構成する細胞の数を推定することが可能である。内細胞塊Ｉが占める領域の取り扱いについても、二次元マップのケースと同様にすることができる。

以上のように、内部の空洞を覆う層の細胞の評価に際し、三次元像から求められる厚さプロファイルにおいて周囲より厚い突起部位を検出し、これを個々の細胞に対応するものとみなすことで、細胞の定量的評価を容易にすることが可能である。対象物の外形が概ね球形である場合、画像データを極座標空間で表すことにより、処理を容易にすることができる。

なお、本発明は上記した実施形態に限定されるものではなく、その趣旨を逸脱しない限りにおいて上述したもの以外に種々の変更を行うことが可能である。例えば、上記実施形態の画像処理装置１は、試料ＳをＯＣＴ撮像および光学顕微鏡撮像する機能、ならびに撮像データから出力画像を作成し出力する機能を有するものである。しかしながら、本発明の画像処理方法は、自身は撮像機能を持たず、撮像機能を有する他の装置での撮像により得られた撮像データを取得したコンピュータ装置によって実行することも可能である。これを可能とするために、図３の各処理ステップのうちステップＳ１０１以外をコンピュータ装置に実行させるためのソフトウェアプログラムとして、本発明が実施されてもよい。

このようなプログラムの配布は、例えばインターネット等の電気通信回線を介してダウンロードする形式によって行うことが可能であり、また当該プログラムを記録したコンピュータ読み取り可能な記録媒体を配布することによっても可能である。また、既存のＯＣＴ撮像装置にインターフェースを介してこのプログラムを読み込ませることで、当該装置により本発明を実施することも可能となる。

また、上記実施形態では三次元画像データを極座標データに変換した上で各種の処理を実行しているが、このような座標変換を行わなくても、上記と同様の突起部位の検出は可能である。例えば厚さ検出にDistance Transform関数を用いる場合、必ずしも胚の中心（重心）を設定する必要はないから、直交座標系で表された栄養外胚葉の三次元画像データから直接厚さを求めることが可能である。そして、図１０（ｂ）に示される三次元マップを用いる場合、内側表面または外側表面の位置が特定されている限り、ＸＹＺ直交座標空間内でマップ作成が可能である。したがって、これらを組み合わせた処理であれば、極座標系への変換は不要であると言える。

また例えば、上記実施形態では、胚Ｅの三次元像から透明帯と栄養外胚葉とを分離した上で、栄養外胚葉における厚さ分布を求めている。しかしながら、例えば透明帯の厚さが一定と見なせるのであれば、透明帯と栄養外胚葉とを分離することなくその厚さ分布を求めても、それに現れる厚さの変化は実質的に栄養外胚葉の厚さを反映したものとなる。

また、上記実施形態は、本発明を胚盤胞期の胚の評価に適用したものである。しかしながら、本発明の適用範囲はこれに限定されず、内部に空洞を有する各種の細胞塊の観察・評価に適用可能である。例えば腸管上皮細胞のオルガノイドは上記と同様の空洞を内部に有し、その表面は細胞により形成された層に覆われている。このような細胞層の評価にも、本発明を好適に適用することができる。

以上、具体的な実施形態を例示して説明してきたように、本発明に係る画像処理方法においては、例えば、細胞塊の内部の一点を原点とする動径に沿った細胞塊の厚さを求めるように構成されてもよい。細胞塊が球体または回転楕円体に近い外形形状を有している場合、このように動径方向に厚さを求めることで処理を容易にすることができる。

また例えば、三次元画像データから細胞塊の各点の位置を極座標空間内の位置として表した極座標データを作成し、極座標データに基づき細胞塊の厚さ分布を求めるように構成されてもよい。特に、細胞塊の重心を極座標空間の原点とすることが好ましい。このような構成によれば、球体または回転楕円体に近い形状の細胞塊についてのデータ処理を容易にすることができる。

この場合において、例えば、極座標空間における動径方向と、当該動径方向において求められた細胞塊の厚さとを関連付けた細胞塊の厚さプロファイルを作成し、これに基づき突起部位を検出するようにしてもよい。このような構成によれば、厚さプロファイルにより示される細胞塊の厚さの変化から、周囲に対し厚さが突出する突起部位を容易に検出することが可能である。

具体的には、例えば厚さプロファイルに対し最大値フィルタ処理を実行し、その前後で値に変化のない位置を突起部位の位置とすることができる。突起部位の位置では細胞塊の厚さが周囲より大きいため、最大値フィルタ処理によっても厚さの値が変化しない。一方、周囲により厚い部分がある領域では、その厚さの影響により、最大値フィルタ処理後の値が変化する。このことを利用して、周囲よりも大きな厚さを有する部位を検出することができる。

例えば、細胞塊の内側表面または外側表面に対応する位置に、当該位置における細胞塊の厚さに応じた輝度値を有する画素を配した三次元マップにより厚さプロファイルを表すことが可能である。また例えば、極座標空間内で細胞塊を球体近似した近似球の表面の各位置に、当該位置を通る動径方向における細胞塊の厚さに応じた輝度値を有する画素を配する三次元マップにより厚さプロファイルを表すことが可能である。これらの場合、注目画素を中心とする球形領域に対する空間フィルタ処理として最大値フィルタ処理を実行することができる。

一方、極座標データが表す２つの偏角を座標軸とする座標平面に、当該２つの偏角により極座標空間内で特定される一の動径方向における細胞塊の厚さに応じた輝度値を有する画素を配した二次元マップにより厚さプロファイルを表すことも可能である。この場合、二次元マップ上の注目画素を中心とする円形領域に対する二次元フィルタ処理として最大値フィルタ処理を実行することができる。なお、三次元空間での厚さ分布を二次元マップにより表現する方法としては、例えば擬円筒図法を用いることが可能である。

また例えば、極座標空間を原点を通る平面で複数に分割し、分割された空間のそれぞれについて個別に二次元マップの作成および前記突起部位の検出を実行してもよい。三次元空間内の物体の形状を二次元マップで表す場合、マップ上の物体に歪みが生じることは避けられない。空間を分割して個別にマップを作成することで、特定の箇所に大きな歪みが生じるのを回避することができ、厚さ分布の検出精度を高めることができる。

本発明に係る画像処理方法は、検出された突起部位の数を計測する工程をさらに備えていてもよい。このように定量的情報を求めることで、ユーザに対しより有用な情報を提供することが可能となる。例えば対象物の細胞塊が胚盤胞期の胚や腸管上皮細胞オルガノイドである場合、その表面を構成する細胞の計数に、本発明を適用することが可能である。

この発明は、内部に空洞を有する細胞塊の観察および評価作業を支援するのに好適なものであり、例えば培養された胚（受精卵）の状態を評価する作業を支援して、生殖補助医療における、より妊娠成功率が高い良好な胚を選択する目的に利用することができる。

１ 画像処理装置
１０ 保持部
２０ 撮像部
３０ 制御ユニット
３３ 信号処理部
３４ 撮像制御部
Ｅ 胚（細胞塊）
Ｓ 試料
Ｓｉ 内側表面
Ｓｏ 外側表面

Claims (14)

1. 内部に空洞を有する細胞塊の画像を処理対象とする画像処理方法であって、
   前記細胞塊を光干渉断層撮像して、前記細胞塊の三次元像を表す三次元画像データを取得する工程と、
   前記三次元画像データに基づき、前記細胞塊のうち外部空間に面する外側表面と前記空洞に面する内側表面との間の前記細胞塊の厚さを位置ごとに表す厚さ分布を取得する工程と、
   前記厚さ分布の取得結果に基づき、周囲の領域に比較して厚さが突出した突起部位を検出する工程と
   を備える画像処理方法。
2. 前記細胞塊の内部の一点を原点とする動径に沿った前記細胞塊の厚さを求める請求項１に記載の画像処理方法。
3. 前記三次元画像データから、前記細胞塊の各点の位置を極座標空間内の位置として表した極座標データを作成し、
   前記極座標データに基づき前記細胞塊の厚さ分布を求める請求項１または２に記載の画像処理方法。
4. 前記細胞塊の重心を前記極座標空間の原点とする請求項３に記載の画像処理方法。
5. 前記極座標空間における動径方向と、当該動径方向において求められた前記細胞塊の厚さとを関連付けた前記細胞塊の厚さプロファイルに基づき前記突起部位を検出する請求項３または４に記載の画像処理方法。
6. 前記厚さプロファイルに対し最大値フィルタ処理を実行し、その前後で値に変化のない位置を前記突起部位の位置とする請求項５に記載の画像処理方法。
7. 前記内側表面または前記外側表面に対応する位置に、当該位置における前記細胞塊の厚さに応じた輝度値を有する画素を配することで前記厚さプロファイルを表す三次元マップを作成し、
   前記最大値フィルタ処理を、注目画素を中心とする球形領域に対する空間フィルタ処理として実行する請求項６に記載の画像処理方法。
8. 前記極座標空間内で前記細胞塊を球体近似した近似球の表面の各位置に、当該位置を通る動径方向における前記細胞塊の厚さに応じた輝度値を有する画素を配することで前記厚さプロファイルを表す三次元マップを作成し、
   前記最大値フィルタ処理を、注目画素を中心とする球形領域に対する空間フィルタ処理として実行する請求項６に記載の画像処理方法。
9. 前記極座標データが表す２つの偏角を座標軸とする座標平面に、当該２つの偏角により前記極座標空間内で特定される一の動径方向における前記細胞塊の厚さに応じた輝度値を有する画素を配することで前記厚さプロファイルを表す二次元マップを作成し、
   前記最大値フィルタ処理を、前記二次元マップ上の注目画素を中心とする円形領域に対する二次元フィルタ処理として実行する請求項６に記載の画像処理方法。
10. 擬円筒図法により前記二次元マップを作成する請求項９に記載の画像処理方法。
11. 前記極座標空間を原点を通る平面で複数に分割し、前記二次元マップの作成および前記突起部位の検出を、分割された空間のそれぞれについて個別に実行する請求項９または１０に記載の画像処理方法。
12. 検出された前記突起部位の数を計測する工程をさらに備える請求項１ないし１１のいずれかに記載の画像処理方法。
13. 内部に空洞を有する細胞塊を光干渉断層撮像した、前記細胞塊の三次元像を表す三次元画像データに基づき、前記細胞塊のうち外部空間に面する外側表面と前記空洞に面する内側表面との間の前記細胞塊の厚さを位置ごとに表す厚さ分布を取得する工程と、
    前記厚さ分布の取得結果に基づき、周囲の領域に比較して厚さが突出した突起部位を検出する工程と
    を、コンピュータに実行させるためのプログラム。
14. 請求項１３に記載のプログラムを記録した、コンピュータ読み取り可能な記録媒体。

Images