JPWO2007069414A1 - 線状の形態を有する細胞等を解析する方法及び神経細胞解析方法並びにそれら方法を実行する装置及びプログラム - Google Patents

Abstract

神経細胞の３次元画像に基づき、神経細胞の形態解析を自動的に行う技術を提案する。まず、スケールスペース法を用いて樹状突起のトレースを行う。この際、σコンボリューションぼかしを利用して凹凸を減らし、樹状突起の中心線を同定する。次に、ヘステンソル法を用いて負の曲率を探す。全ての座標軸で負の曲率となる部分を「頭部」が占める領域であると判断する。この領域（頭部）の中心をスパインの位置とする。スパインの頭部は、楕円体と見なして近似を行う。これによって、楕円体の短径、中径及び長径をそれぞれ算出する。また、スパインの位置からそのスパインに最も近い樹状突起に垂線を降ろすことによって、その垂線を柱部分と認定する。このようにして得られた樹状突起と、スパイン頭部及び柱部とを合わせて最終的な神経細胞の形態が得られる。

Description

本発明は、線状の形態を有する細胞や、組織、器官を解析する方法に関する。特に、本発明は、神経細胞の解析方法に関する。特に、神経細胞の樹状突起の分岐の様子や、スパインの位置や形態、等を特定し、解析する方法に関する。

また、これらの方法を実行する装置、及び、その方法を実行するプログラムに関する。

従来から、神経細胞の解析は大きな労力を必要としていた。例えば、顕微鏡画像を基にした神経細胞の可視化解析は、熟練研究者の目で１個１個確認していく手作業によって、単一の神経細胞の形態情報をコンピュータに入力する等の方法がもっぱら採用されている。このような作業は、可視化に際して、可視化のための簡単なソフトウェアを用いる場合も多いが、結局、人手による作業が多くを占めているので、その作業量は膨大なものとなりがちであった。

例えば、１枚の画像解析を完了するのに数日を要する場合も珍しくはなかった。

これは、１個の神経細胞には、数十程度の樹状突起が存在し、各樹状突起上に数万程度のスパイン（シナプス後部）が存在するという複雑な形態を神経細胞がそなえているためである。このような複雑な構造のため、従来の画像処理アルゴリズムでは、自動的に神経細胞の形態を解析することは困難である。

先行特許文献の例
例えば、下記特許文献１には、線状構造を有する細胞体の評価方法が記載されている。この特許文献１に記載の方法は、画像データから、神経突起などの伸長度合いを取得する方法である。この方法は、細胞体の個数及び面積と、線状構造の長さ及び面積と、を画像処理によって求めることを特徴とするものである。

また、下記特許文献２には、蛍光標識を有するレポータ分子を用いて神経突起の成長に特異的に影響する化合物をスクリーニングする方法が、記載されている。この方法に用いられる光学分析システムや、試薬等も同特許文献２に記載されている。

また、下記特許文献３には、細胞スクリーニングの方法が記載されている。この方法においては、顕微鏡と、その画像をデジタル画像化するシステムと、細胞物質を区分するためのソフトウェア等が用いられている。

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

このスパインは、神経の活動に伴い増減し、「記憶」において重要な役割を果たすことが本願発明者らの研究により明らかになっている。このスパインを観察・解析することは、記憶の動作を解明し、種々の脳疾患の治療法の開発に寄与するものであり、医学上、重要な意味を持っている。

さて、このスパインは、１ニューロンあたり、数万個存在するので、それを人力で解析するのは極めて過大な労力が必要であった。特に、単にスパインの個数だけでなく、その形態まで観察し解析するためには熟練の技術者が何日も作業を続ける必要があった。

したがって、神経細胞の解析、特にスパインの増減まで含めた神経細胞の形態を高速に解析する手法が従来から望まれていた。

本願発明は、係る社会的な要請に鑑みなされたものであり、その目的は、神経細胞の３次元画像に基づき、神経細胞の形態解析を自動的に行う技術を提案することである。また、この技術を応用し、線状の形態を有する細胞等の形態を解析する方法・装置等を提供することである。

Ａ．方法
（１）本発明は、上記課題を解決するために、線状の形態を有する細胞又は組織又は器官である解析対象物の３次元画像に基づき、前記対象物を解析する方法において、前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかしステップと、前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度ステップと、前記第１の中心度に基づき、前記対象物の線状の形態に近似した図形を求める線状形態フィッティングステップと、を含むことを特徴とする線状の形態を有する細胞又は組織又は器官を解析する方法である。

このような構成によって、線状の形態を有する細胞等の形態を解析することができる。例えば、神経細胞や、（毛細）血管等の解析を好適に行うことができる。

（２）神経細胞の３次元画像に基づき、神経細胞の形態を解析する方法において、前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかしステップと、前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度ステップと、前記第１の中心度に基づき、神経細胞の樹状突起に近似した図形を求める樹状突起フィッティングステップと、を含むことを特徴とする神経細胞解析方法である。

このような構成によって、神経細胞の形態、特に樹状突起の形態を自動的に解析することができる。

（３）また、本発明は、上記（２）記載の神経細胞解析方法において、前記第１のぼかしステップは、前記第１のスケールσの値を変化させてスケールスペース変換を行い、各スケールσに対する複数のぼかし画像群を生成することを特徴とする神経細胞解析方法である。

（４）また、本発明は、上記（３）記載の神経細胞解析方法において、前記第１のぼかしステップは、前記第１のスケールσの値を０．２μｍ−５．０μｍの範囲で変化させてスケールスペース変換を行い、各スケールσに対する複数のぼかし画像群を生成することを特徴とする神経細胞解析方法である。

樹状突起の形態を把握するために、σの値として０．２μｍ−５．０μｍの範囲が好適である。

（５）また、本発明は、上記（２）記載の神経細胞解析方法において、前記第１の中心度ステップは、境界面の存在度を中心度として求めることを特徴とする神経細胞解析方法である。

中心度として境界面の存在度を用いたので、対象物の表面を良好に把握することができる。つまり、ここで言う境界面とは、対象物（樹状突起）の表面を意味する。

（６）また、本発明は、上記（２）２記載の神経細胞解析方法において、前記第１の中心度ステップは、はぎ取り法を用いて中心線を求めることを特徴とする神経細胞解析方法である。

従来のはぎ取り法を用いても、精度は多少落ちるものの、樹状突起の形態を同様に求めることが可能である。

（７）また、本発明は、上記（２）記載の神経細胞解析方法において、前記樹状突起フィッティングステップは、前記求めた中心度に基づき、スネークを樹状突起にフィッティングさせることを特徴とする神経細胞解析方法である。

Ｂ１．装置
（８）また、本発明は、線状の形態を有する細胞又は組織又は器官である解析対象物の３次元画像に基づき、前記対象物を解析する装置において、前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかし手段と、前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度手段と、前記第１の中心度に基づき、前記対象物の線状の形態に近似した図形を求める線状形態フィッティング手段と、を含むことを特徴とする線状の形態を有する細胞又は組織又は器官を解析する装置である。

Ｃ１．プログラム
（９）また、本発明は、コンピュータを、線状の形態を有する細胞又は組織又は器官である解析対象物の３次元画像に基づき、前記対象物を解析する装置として動作させるプログラムにおいて、前記コンピュータに、前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかし手順と、前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度手順と、前記第１の中心度に基づき、前記対象物の線状の形態に近似した図形を求める線状形態フィッティング手順と、を実行させることを特徴とするプログラムである。

Ａ２．方法
（１０）また、本発明は、上記（１）記載の線状の形態を有する細胞・組織又は器官を解析する方法において、前記３次元画像に対して、所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかしステップと、前記第２のぼかしステップによるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度ステップと、前記第２の中心度に基づき、前記線状の形態の上に存在する棘状部位の検出を行う棘状部位検出ステップと、を含むことを特徴とする線状の形態を有する細胞・組織又は器官を解析する方法である。

ここで、棘状部位とは、棘状であればどのようなものでも含まれる。例えば、神経細胞上のスパインが典型的な一例である。

（１１）また、本発明は、上記（２）記載の神経細胞解析方法において、前記３次元画像に対して、前記第１のスケールσと異なる所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかしステップと、前記第２のぼかしステップによるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度ステップと、前記第２の中心度に基づき、神経細胞のスパイン頭部の中心位置を検出するスパイン頭部検出ステップと、を含むことを特徴とする神経細胞解析方法である。

（１２）また、本発明は、上記（１１）記載の神経細胞解析方法において、前記第２のスケールσの値を０．１μｍ−０．３μｍの範囲で変化させてスケールスペース変換を行い、各スケールσに対する複数のぼかし画像群を生成することを特徴とする神経細胞解析方法である。

（１３）また、本発明は、上記（１１）記載の神経細胞解析方法において、前記スパイン頭部検出ステップは、前記第２の中心度の絶対値がピークを取る位置を、スパインの頭部の中心位置と判断することを特徴とする神経細胞解析方法である。

（１４）また、本発明は、上記（１１）記載の神経細胞解析方法において、前記スパイン頭部検出ステップは、前記樹状突起から所定の距離以内の位置を、スパインの頭部の中心位置と判断することを特徴とする神経細胞解析方法である。

（１５）また、本発明は、上記（１１）記載の神経細胞解析方法において、前記求めたスパイン頭部の中心位置から、最も近い前記樹状突起へ、垂線を降ろし、その垂線をスパインの柱部として求めるスパイン柱部検出ステップ、を含むことを特徴とする神経細胞解析方法である。

（１６）また、本発明は、上記（１１）記載の神経細胞解析方法において、前記スパイン検出ステップは、ヘステンソルの固有値が全て負である領域をスパインの頭部の存在領域として求めることを特徴とする神経細胞解析方法である。

（１７）また、本発明は、上記（１６）記載の神経細胞解析方法において、前記スパイン検出ステップは、前記求めた存在領域を、楕円体で近似し、スパイン頭部の近似図形を求めることを特徴とする神経細胞解析方法である。

Ｂ２．装置
（１８）また、本発明は、上記（８）記載の線状の形態を有する細胞・組織又は器官を解析する装置において、前記３次元画像に対して、所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかし手段と、前記第２のぼかし手段によるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度手段と、前記第２の中心度に基づき、前記線状の形態の上に存在する棘状部位の検出を行う棘状部位検出手段と、を含むことを特徴とする線状の形態を有する細胞・組織又は器官を解析する装置である。

Ｃ２．プログラム
（１９）また、本発明は、上記（９）記載の、コンピュータを、線状の形態を有する細胞・組織又は器官を解析する装置として動作させるプログラムにおいて、前記３次元画像に対して、所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかし手順と、前記第２のぼかし手順によるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度手順と、前記第２の中心度に基づき、前記線状の形態の上に存在する棘状部位の検出を行う棘状部位検出手順と、を実行させることを特徴とするプログラムである。

D. 神経細胞解析方法の他の例
（２０）また、本発明は、上記（２）記載の神経細胞解析方法において、前記樹状突起フィッティングステップは、前記求めた中心度に基づき、中心度が最大の点を選択する第１選択ステップと、前記中心度が最大の点から所定距離離れた点の中で、前記中心度が最大の点を選択する第２選択ステップと、前記選択した２個の点を結び、初期スネークを構成するステップと、前記構成した初期スネークの端点から所定距離離れた点の中から、最大の中心度を持つ点を選択し、既存の端点と結び、前記スネークを成長させる成長ステップと、を含むことを特徴とする神経細胞解析方法である。

（２１）また、本発明は、上記（２０）記載の神経細胞解析方法において、前記成長ステップにおいて、選択する点は、前記既存の端点との距離が所定の距離にある点であることを特徴とする神経細胞解析方法である。

（２２）また、本発明は、上記（２０）又は（２１）記載の神経細胞解析方法において、前記成長ステップにおける選択する点を見つけられる場合、前記成長ステップにおける成長を続行し、前記選択する点を見つけられない場合は、前記成長ステップにおける成長を完了する制御ステップ、を含むことを特徴とする神経細胞解析方法である。

（２３）また、本発明は、上記（２２）記載の神経細胞解析方法において、前記成長ステップにおける選択する点を見つけられる場合とは、前記選択する点を、画像中において既にスネークの内部と認定されている領域を除いた領域から見つけられる場合であることを特徴とする神経細胞解析方法である。

（２４）また、本発明は、上記（２０）、（２１）、（２２）又は（２３）記載の神経細胞解析方法において、前記成長させたスネークの各ノードで直径を取得し、スネーク内部の領域を確定させる確定ステップ、を含むことを特徴とする神経細胞解析方法である。

（２５）また、本発明は、上記（２０）、（２１）、（２２）、（２３）又は（２４）記載の神経細胞解析方法において、前記成長させたスネークが、（ａ）３個以上のノードを有すること、（ｂ）全ノードの平均中心度が、所定のそれまでに確定した樹状突起の平均中心度の所定の比率以上であること、の２条件をともに満たす場合に樹状突起として確定し、いずれかの条件を満たさない場合には、そのスネークを破棄し、前記第１選択ステップに戻り、中心度が最大の点を選択する動作を再び開始する繰り返し制御ステップ、を含むことを特徴とする神経細胞解析方法である。

以上述べたように、本発明によれば、線状の形態を有する細胞等の解析を効率的に実行することができる。特に、対象が神経細胞の場合は、その樹状突起の様子とともに、樹状突起上に表れるスパインの形態も効率的に解析することができる。

神経細胞の形状の説明図である。 ３次元画像データの説明図である。 種々のスパインの形態の説明図である。 ストレス・ステロイド作用に伴うスパインの数の変化（密度の変化）の説明図である。 従来の画像解析の様子を説明する説明図である。 本実施の形態における特徴的な３ステップの概念図である。 σコンボリューションによるノイズ除去の様子を示す説明図である。 中心度の画像である。 シリンダ図形を連結させて樹状突起を表現した例を示す説明図である。 はぎ取り法の説明図である。 ヘステンソルの固有値の積の説明図である。 スパイン頭部の候補を選ぶ際に、擬陽性の候補を除去する様子を示す説明図である。 主軸の説明図である。 勾配ベクトルの説明図である。 本実施例の動作を表すフローチャートである。 実施例2の動作を表すフローチャートである。

符号の説明

１０ ゴミ
１２ 擬陽性

以下、本発明の好適な実施の形態を図面に基づき説明する。

１．神経の形態
図１に、神経細胞の形態の説明図が示されている。本実施の形態においては、神経細胞の形態は、対象となる神経細胞に蛍光色素をインジェクションし、これを共焦点顕微鏡によって観察することによって、神経細胞の画像データが得ている。この画像データは、１０２４×１０２４ドットの画像であり、１ドットが０．０４３μｍの大きさである。

そして、共焦点顕微鏡で観察する（撮影する）平面の位置を、奥行き方向（平面と垂直方向）に変化させることによって、奥行き方向に０．４５μｍずつ位置が異なる平面画像データを得ることができる。

この様子が、図２に示されている。本実施の形態では、０．４５μｍずつずれた位置の平面画像データを例えば３０枚撮影している。これら３０枚の１０２４×１０２４ドットの画像データ群は、神経細胞の形態を３次元的に表現しているので、これら３０枚の画像データ群を本実施の形態では、「３次元画像データ」と呼ぶ。

さて、図１に示すように、神経細胞からは、周囲に向かって樹状突起が伸展しており、複雑な形態をなしている。そして、この樹状突起上には、スパイン（棘突起と呼ぶ場合がある）と呼ばれる突起が発生することが知られている。

スパインの形態を示す図が図１の（２）に示されている。このように、スパインは、数ミクロン程度の小さな突起であり、単一の神経に対しておよそ１万個程度あることが知られている。また、密度は、１スパイン／１μｍ程度である。

このスパインは、本願発明者らの研究によって、記憶に関して重要な働きをしていることが判明しており、その発生の様子を観察することは、記憶のメカニズムを探る上で極めて重要な事項となっている。ところが、このスパインは、１個の神経細胞あたり１万個程度存在し、スパインの形態を画像データから人手で解析するのは極めて繁雑な作業が必要であった。

本実施の形態では、神経細胞の形態を効率的に求めることができる手法を説明する。特に、樹状突起の形態、スパインの形態、を含む神経細胞の形態を効率的に求める手法を説明する。

２．スパインの形態
図３には、種々のスパインの形態を説明する図が示されている。この図に示すように、スパインは、ｔｈｉｎ、ｓｔｕｂｂｙ、ｍｕｓｈｒｏｏｍ、ｆｉｌｏｐｏｄｉｕｍ、の４種類の形態が知られている。

図４には、ストレス・ステロイド作用に伴うスパインの数の変化（密度の変化）の説明図が示されている。この図のグラフにおいて、黒はｍｕｓｈｒｏｏｍ、斜線はｔｈｉｎ、白はｆｉｌｏｐｏｄｉｕｍ、縦縞はｓｔｕｂｂｙの密度の変化をそれぞれ表す。

このように、ストレスを加えることによって、スパインの数（密度）が変化することが知られている。

３．従来のスパイン解析
現在、神経細胞の形態を調べる標準的な手法は、神経細胞の画像データから、手作業によって神経細胞の形態、特にスパインの個数等を求める手法である。典型的には、Ｎｅｕｒｏｌｕｃｉｄａ（商品名）という名のプログラムを用いて神経細胞の画像解析を行う。 この画像解析の様子が図５に示されている。図５（１）は、共焦点顕微鏡で得られた原画像である。この原画像に対して、操作者が目視でスパインの位置を決定し、スパインの柱部の長さと太さを決定していた。この様子が図５（２）の写真に示されている。

しかしながら、スパインは、簡略化すれば、概ね、頭部（ヘッド）と、頭部と樹状突起を結ぶ柱部と、から構成されていることが知られている（図５（３）参照）。

望ましくはこの図５（３）に示すように、楕円体である頭部の幅ｗ１、ｗ２と、同じく頭部の長さ（楕円体の長辺）ｔと、柱部の長さｈと、この柱部と樹状突起とのなす角度θｌと、を各スパイン毎に求めることが好ましい。

上記従来の手法では、長さｔと、幅ｗ（ｗ１とｗ２の平均値と考えられる値）しか判定できなかった。また、頭部の楕円体のパラメータ（従来の場合はｗ、ｔ）を得るまで、計算をやり直す必要があり、操作者の操作労力が過大なものとなりがちであった。また、手作業の割合が多いため、膨大なデータには対応が困難であるという問題点もあった。

本実施の形態では、これらのパラメータの内、ｗ１、ｗ２、ｔ、ｈを自動的に算出する例を示す。

４．本実施の形態における基本アルゴリズムとその特徴
本実施の形態においては、共焦点顕微鏡にて撮影したの画像（例えば１０２４ドット×１０２４ドット）を複数画面（例えば、２０面−３０面）用意し、３次元画像データを構成している。

これら３次元画像データを入力し、神経細胞の形態を迅速に得ることができる方法を提案する。本実施の形態における方法は、以下の３個の特徴的なステップを含むものである。この３ステップの概念図が図６に示されている。

（ステップ１）
樹状突起のトレース：スケールスペース法を用いて樹状突起のトレースを行う。この際、σコンボリューションぼかしを利用して凹凸を減らし、樹状突起の中心線を同定する。このステップ１の概念図が図６のステップ（１）に示されている。

（ステップ２）
スパイン位置の検出：ヘステンソル法を用いて負の曲率を探す。負の曲率となる部分が、ラグビーボール状の「頭部」に相当すると考えられるので、全ての座標軸で負の曲率となる部分を「頭部」であると判断する。本ステップ（２）におけるスパイン位置とは、この「頭部」の位置である。このステップ２の概念図が図６のステップ（２）に示されている。

（ステップ３）
スパインの頭部は、楕円体（三軸不等楕円体）で近似できると考えられるので、楕円体と見なして近似を行う。これによって、短径、中径及び長径をそれぞれ算出する。なお、主軸方向の径の内、最も長い径を長径（ｔ）と呼び、２番目に長い径を中径（ｗ２）、最も短い径を短径（ｗ１）と呼ぶ。図５の（３）を参照されたい。

また、スパインの柱部は、スパインからそのスパインに最も近い樹状突起に垂線を降ろすことによって、その垂線を柱部と認定する。このステップ３の概念図が図６のステップ（３）に示されている。

このようにして得られた樹状突起と、スパイン頭部及び柱部とを合わせて最終的な神経細胞の形態が得られる。

以下、本実施の形態の神経細胞解析方法の動作を、上記各特徴的なステップの動作原理を中心に説明する。

４−１． ステップ１：樹状突起のトレース
４−１−１．σコンボリューション
本実施の形態においては、スケール・スペースを用いて尺度非依存的な解析法を採用している。すなわち、対象物の大きさにあったスケールσを選択している。例えば、樹状突起をトレースする場合には、スケールσ＝０．２−５．０μｍの範囲で設定する。また、もしスパインが対象である場合には、σ＝０．１−０．３μｍの範囲で設定する。この結果、このスケールσより小さい凸凹構造を「σコンボリューション」でぼかしてノイズを除去するのである。

本ステップ１では、樹状突起のトレースを行うので、スケールσ＝０．２−５．０μｍに設定している。

このσコンボリューションの式は、以下の通りである。

ここで、Ｉは、画像を構成する各画素の輝度であり、Ｋ（σ）とのコンボリューションを取ることによって、ぼかしを与えている。Ｋ（σ）の式は上記式の通りである。ここで、Ｎは、正規化のための定数であって、カーネル関数Ｋ（σ）を空間全体で積分した結果が１となる用に定めた数である。例えば、２次元の場合は、Ｎ＝１／（２πσ^２）である。

なお、上記数１では、ｘ軸方向のみが記されているが、同様の処理をｙ軸、ｚ軸についても行う。

このσコンボリューションによるノイズ除去の様子が図７に示されている。図７（１）は、原画像データであり、図７（２）は、σ＝０．５μｍに設定した場合のσコンボリューションの結果である。図７（３）は、σ＝１．０μｍに設定した場合のσコンボリューションの結果である。図７（４）は、σ＝２．０μｍに設定した場合のσコンボリューションの結果である。本実施の形態においては、まず、神経細胞の樹状突起をトレースするためにノイズ除去をするので、σ＝０．２μｍ−５．０μｍの範囲でσを設定している。

実際には、σの値は、０．２μｍ−５．０μｍの範囲でステップ的に変化させて複数のぼかし画像を作成する。例えば、σ＝０．５μｍ、σ＝０．６μｍ、σ＝０．７μｍ、・・・のようにσを変化させて各σ毎にぼかし画像を生成する。後述する実施例では、より詳細なぼかし画像群の生成動作を説明している。

なお、ぼかし画像群は２種類作成している。第１に、上述した樹状突起をトレースするためのぼかし画像群である。第２に、スパイン頭部を検出するためのぼかし画像群である。

このスパイン頭部を検出するための第２のぼかし画像群は、変化させるσの値の範囲が０．１−０．３μｍであることのみが異なり、その他は樹状突起用の第１のぼかし画像群と同様にして作成する。

４−１−２． 中心度変換
上記スケールスペース変換（σコンボリューション）でノイズ除去をした後、中心度変換を行って樹状突起のトレース、及びスパイン頭部の検出を行う。

樹状突起のトレースの場合は、境界面の存在度を中心度として用いている。一方、スパイン頭部の検出の場合は、中心度としてヘステンソルの固有値の積を用いている。以下、樹状突起のトレースについて説明する。スパイン頭部の中心度の算出は、次節４−２において説明する。

樹状突起のトレース
まず、スケールスペース変換後のぼかし画像群の各画像に対して、各画素のヘステンソルを求め、その固有ベクトルを求める。この固有ベクトルの内、最小の固有値（負の値）に属するベクトルの方向（両側）に上述したσだけ「進んだ」位置（画素）の輝度の第１の勾配ベクトルを求める。同様にσだけ「後退した」位置（画素）の輝度の第２の勾配ベクトルを求める。

そして、上記第１の勾配ベクトルと、上記固有ベクトルを正規化したベクトルと、の内積を求める。この値に−１を掛けたものを「第１のスケールσでの境界度」と呼ぶ。同様に、上記第２の勾配ベクトルと、上記固有ベクトルを正規化したベクトルと、の内積を求める。この値を「第２のスケールσでの境界度」と呼ぶ。

そして、「第１のスケールσでの境界度」と「第２のスケールσでの境界度」とのいずれか小さい方（正の値）をその画素における「中心度」とする。

この中心度は、境界面の存在度を表す尺度であり、境界面とは、対象物すなわち樹状突起の表面を言う。

さて、この中心度の輝度（値）が極大となる位置をトレースしていけば、樹状突起の幹線（中心線）と幅（太さ）をトレースすることができる。なお、幅（太さ）は、極大中心度を与えるスケールσの値である。

この中心度の画像が図８に示されている。図８（１）は、原画像データであり、図８（２）は、σ＝０．５μｍの場合に、各画素について求めた中心度の画像（中心度画像と呼ぶ）である。同様に、図８（２）は、σ＝１．５μｍの場合に、各画素について求めた中心度画像である。同様に、図８（３）は、σ＝１．５μｍの場合に、各画素について求めた中心度画像である。図８（４）は、σ＝２．０μｍの場合に、各画素について求めた中心度画像である。

このように、複数のぼかし画像群のそれぞれに対して、中心度画像をそれぞれ求める。

これらの図において、中心度が大きい値の画素（明るい画素）は、半径σのシリンダ図形がその画素を通っていることを示している。

すなわち、図８（２）−図８（４）の画像中、明るさの極大値をトレースした線が、樹状突起の中心線である。そして、その中心線を中心として半径σのシリンダ図形がその画素を通っていることを示している。

さて、シリンダ図形とは、簡単に言えば、円錐を所定の断面で切断した図形モデルであるが、詳しくは後述する。本実施の形態では、樹状突起をこのシリンダ図形を連結させることによって表現（近似）している。実際にシリンダ図形を連結させて樹状突起を表現した例を図９に示す。

本実施の形態においては、エネルギー関数として、中心度（medialness）と歪みとの和から成るエネルギー関数を定義し、このエネルギー最小問題を解くことによって実際の樹状突起にフィットさせている。

なお、フィット（又はフィッティング）とは、簡単に言えば近似することを言う。より詳細に言えば、「現実に存在する図形に、理想的な基本図形を対応させること」を意味する。また、「誤差を最小とするような、対応する基本図形を与える」ことでもある。

さて、本実施の形態におけるエネルギー関数の式を以下に示す。

本実施の形態では、このエネルギー関数を、スネーク関数と呼んでいる。

スネークとは、位置、太さ、接続情報からなる節点（スナクセル）の集合である。ここで、接続情報とは、ある節点が他のどの節点とつながっているかを示す情報である。位置と太さを空間的パラメータと呼び、接続情報を位相的パラメータと呼ぶ。本実施の形態においては、接続情報は、近接する高々二つの節点と自動的に接続することで作成するので、フィッティングは空間的パラメータについてのみ行う。

空間的パラメータは位置（三次元）と太さ（一次元）の複合体なので、四次元の空間で考えるのが妥当である。そして、この四次元空間がスケールスペース（尺度空間）である。中心度を求めることは、この四次元空間内に「線図形の存在し得る度合い」をプロットすることに相当する。

上記式中、Ｅｅｘｔ（ｉ） は、ｉ番目の節点の位置と太さにおける中心度の符号を反転させたものである。この値が小さいということは、その節点のパラメータが指定する位置に、指定された太さの直線図形が通っていることを意味する。次に、Ｅｉｎｔ（ｉ）は、ｉ番目の節点の内部エネルギーを表す。これは、スネークを必要以上に曲げないようにするためのエネルギーで、この節点におけるスネークの曲率を用いる。この値が小さいということは、その節点においてスネークが直線的であることを意味する。αは分配定数で、フィッティングの最中には変化しない。αが大きい場合には中心度エネルギーの方を重視し、αが小さい場合には内部エネルギーの方を重視してフィッティングを行う。

Ｅｅｘｔ（ｉ）は、上で述べたように、ｉ番目の節点の位置、太さにおける中心度（但し符号は反転）である。スネークを必要以上に曲げないという拘束条件下で、この符号反転した中心度（エネルギー）を最小化することで、スケールスペース中の中心度の高い部分を通過するスネークを求めることができる。これがスネークのフィッティングである。フィッティングにおける変化パラメータは、空間的パラメータ（節点の位置と太さ）である。

中心度変換の変形例
樹状突起上のスパインが比較的大きい場合に、上で述べた中心度変換を行うと、スパインも棘として残ってしまうのでトレースが困難になる場合がある。そこで、従来から使用されているいわゆる「はぎ取り法」を用いて中心線を得ることも好ましい。

はぎ取り法の説明図が図１０に示されている。

はぎ取り法は、物体の中心線を求める手法の一つであり、物体を構成する点（画素）を外側から１画素ずつ削り取って中心線を残す方法である。このはぎ取り法を適用した例が図１０に示されている。図１０（１）は原画像データである。この原画像データを２値化し、その後、外周部から１画素ずつ画素を削り取っていく。その結果が図１０（２）に示されている。

このはぎ取り法は、以下のような欠点が知られている。

（１）境界の凹凸が棘として残り、自動トレースが失敗する場合がある。

（２）はぎ取り結果が、スケールに依存してしまう。

しかし、本実施の形態では、スケールスペース法を用いてσコンボリューションによるぼかしを行っているので、中心度変換にこのはぎ取り法を用いても良好なトレースができる可能性が高い。

はぎ取り法による中心度変換では、得られた中心線に相当する画素に、例えば最大輝度の「２５５」等が代入され、それ以外の画素には最小値である「０」が代入される。これらの値は、正負値を用いて、「１２７」と「−１２８」とを採用しても良い。

４−２．ステップ２：スパイン位置の検出
４−２−１． スパイン頭部の位置
次にスケールσ＝０．２μｍで対角化したヘス・テンソルが３次元で負の固有値を示す位置を求めることによって、球状の部分とその中心を検出する。この球状の部分を検出することは、スパイン頭部の候補を抽出することに他ならない。

なお、スケールσ＝０．２μｍで対角化とは、スケールσ＝ ０．２μｍでぼかし変換を行った後に、求めたヘステンソルを対角化することを意味する。

さて、画像とは、空間の各点に輝度を対応付ける写像である。ヘステンソル（Hessian tensor）は、この輝度分布の二階微分であり、局所的な曲率を表す。

したがって、３次元画像中の負の曲率が、球状に盛り上がったスパインに対応する。本実施の形態では、このことに着目し、ヘステンソルを用いてスパインの存在領域を検出している。

ところで、一般に、ディジタル信号の微分は、差分に置き換えて計算する。例えば、分点公式などを使うのが一般的である。本実施の形態においてもそのような従来の手法を採用することも好ましい。

但し、本実施の形態においては、畳み込み積分の形で書けるぼかし変換を行っているので、画像への微分演算が、畳み込みカーネルへの微分演算に帰着する。すなわち、最初に微分したカーネル関数を用意し、これを畳み込み積分することで、ぼかし変換と微分を同時に行っている。この結果、演算量を節約でき、神経細胞の形態を迅速に求めることが可能となる。

ヘステンソルの式は、以下の通りである。

このヘステンソルを、スケールσ＝０．２μｍで対角化し、固有値を求める。対角後の式を以下に示す。

このようにして求めた固有値λ１、λ２、λ３から、その部分（画素）の曲率が判明する。球状図形の場合は、全ての固有値λ１、λ２、λ３が負であるので、その条件に合致するものだけを取り出して新たな画像を生成する。

すなわち、全ての固有値が負であるものは、その固有値の絶対値（−λ１・λ２・λ３）を画素の値として代入し、それ以外の画素に０を代入した新たな画像を生成する。この画像は、スパイン位置を検出するための中心度画像である。

中心度が０の画素は、そこがスパインの領域ではないことを意味し、中心度が非０の画素は、そこがスパインの存在する領域であることを意味する。

本実施の形態において特徴的なことは、スパイン位置の検出（存在する領域の検出と、スパイン頭部の中心位置の検出）のために、ヘステンソルの固有値の積を中心度として用いていることである。

このようにして、３次元で負の固有値を示す位置（画素）を求めることによって、スパインの頭部に相当する球状の部分（頭部が存在する領域）を求めることができる。また、固有値の積の絶対値のピーク位置からそのスパイン頭部の中心位置を検出することができる。

本実施の形態において特徴的なことは、スパイン頭部の検出において用いた「中心度変換」であり、特に、ヘステンソルの固有値の積を「中心度」として用いていることが大きな特徴である。

なお、スパイン頭部の「中心位置」は、上で述べたように、ピーク値すなわち固有値の積λ１・λ２・λ３の絶対値の極大値を与える点（画素）である。この様子が図１１に示されている。図１１（１）は原画像データであり、図１１（２）は、固有値の積λ１・λ２・λ３の値を示す画像データが示されている。また、図１１（２）においては、固有値の積λ１・λ２・λ３の極大値を取る位置が矢印で示されている。

さて、このようにして得られた球状の部分（固有値の値が全て負の画素）がスパインの頭部の候補となる。３次元画像中の負の曲率が球状のスパインの頭部に由来すると考えられるからである。

ところで、本実施の形態においては、候補を選択する際に、演算の適応範囲を先にフィットした樹状突起の周辺に限定している。このような限定を行うことで、擬陽性１２の候補を除去している。この様子が図１２に示されている。図１２に示すように、本実施の形態では、候補の中心位置（上述したピーク値を与える点（画素））が樹状突起から距離８μｍ以上離れたものは、ゴミ１０と見なしてスパイン頭部の候補から除外している。

なお、ヘステンソルの固有ベクトルは球状の部分を楕円体と考えた場合の主軸を示す。

４−２−２． スパイン頭部及び柱部の詳細な検出
次に、スパインの形態を詳細に解析するため、頭部を楕円体で近似して、その大きさのパラメータを取得する。スパイン頭部に対して、主軸の向きと各軸まわりの慣性モーメント解析から最も妥当な楕円体を得る。なお、主軸は、検出されたスパイン頭部の輝度中心における固有ベクトルの方向である。

このようにして得られた楕円体のパラメータによって、図５で説明したようにｗ１、ｗ２、ｔが求められる。

なお、「主軸」とは、「二次曲面の主軸」のことであり、直感的には、「図形の最も特徴的な方向を指す直交軸」である。例えば、ラグビーボールのような図形の場合、一番長い方向と、それに直交する平面を張る二つの方向が主軸（合計３本）である。

二次元図形の場合には主軸は二本存在するが、三次元図形（実際のスパインの場合）では三本存在する。ちなみに、楕円体も二次曲面の一つである。

主軸の説明図が図１３に示されている。この図においては、一方の矢印が楕円の長軸、他方が短軸を表している。なお、図１３は平面であるので、便宜上、２本の主軸しか示していないが、実際のスパインでは主軸は３本存在する。ヘステンソルを対角化する基底は、この主軸と一致する。

また、主軸の方向に存在する境界面の位置を調べ、その条件で楕円体でフィットする手法も好適である。境界面とは対象物の表面を言い、ここではスパインの頭部の表面を言う。

境界面の位置を調べるとは、元図形に存在した境界面の位置を推定するため、画像（元図形が光学的ボケなどによって劣化したもの）の一階微分と、任意の方向ベクトルｏとの内積を取り、この値が大きい所を「ｏに直交する境界面が存在する場所」とすることである。

この説明図が図１４に示されている。この図１４において、画像の一階微分−ΔＩと、方向ベクトル（ｏ１、ｏ２）が示されている。また、両者の内積の大きさ（ｂ１、ｂ２）が記載されている。この図１４の場合は、ｂ１の大きさの方がｂ２より大きいので、実線のような境界が存在することが判明する。

なお、本特許では、楕円体とは、もちろん回転楕円体だけでなく、三軸不等楕円体も含む意である（図５参照）。

また、スパインの根元、すなわち、柱部が樹状突起と接続する部位は、蛍光が弱く特定することが困難な場合も多い。したがって、本実施の形態では、上述したように、スパインの頭部（中心位置）と樹状突起との最短距離を求めて、これをスパインの柱部と見なしている。このようにして得られた最短距離は、柱部の長さ、すなわちｈとなる（図５参照）。

４−３．
以上のようにして求めたスパイン頭部とその柱部、及び、樹状突起、を合成し、最終的な神経細胞の形態が得られる。

本実施の形態によれば、神経細胞の形態をスパインまで含めた形で得ることができるので、神経活動の研究、疾病の診断等の効率化に大きく資するものである。
「実施例１」

５．実施例１：神経細胞形態自動高速解析プログラム
以上述べたような動作を、本願発明者らは、コンピュータ上のプログラムで実現し、その処理能力を実際に調べた。

処理対象としては、１個の海馬をスライスしたものに、適宜薬物刺激を与えたものを用いた。このスライスを共焦点顕微鏡を用いて撮影した。その際、蛍光色素をインジェクションしておいた。

撮影の解像度は、１０２４×１０２４ドットで、画素あたりのビット数は８ビットであり、２５６階調の輝度を記録した画像データを得た。画素の間隔（画素の大きさ）は、０．０４３μｍである。このような画像データを深さ方向（ｚ方向）に０．４５μｍずつ変化させて３０枚撮影した（図２参照）。

このようにして、１０２４×１０２４ドット（１画素８ビット）の画像データを深さ方向に３０枚得た。これら画像群が、「３次元画像データ」である。

この３次元画像データ（３０枚の画像データ）を、コンピュータに供給し、この３次元画像データに対して、本実施の形態で本願発明者らが作成したプログラムを適用し、神経の形態を出力させた。

本実施例におけるコンピュータの動作（プログラムの動作）を表すフローチャートが図１５に示されている。

５−１．ｚ軸ＭＩＰ画像の作成
まず、図１５のステップＳ１５−１に示すように、上記入力した３次元画像データに基づき、ｚ軸ＭＩＰ画像を作成する。ここで、ｚ軸ＭＩＰ画像とは、所与の断層画像をｚ軸方向に投影してできる二次元画像で、その各点の輝度が、投影線上に並ぶ画素の最大輝度となるような画像のことである。ＭＩＰは、Ｍａｘ Ｉｎｔｅｎｓｉｔｙ Ｐｒｏｊｅｃｔｉｏｎ（最大輝度投影）の略である。

ｚ軸ＭＩＰ画像の概念図が図１６に示されている。この図に示すように、３０枚の画像データの同一の位置にある画素を取り出し（すなわち、３０枚分の３０画素）、その中から最も輝度の高い値を取り出し、その値を、ｚ軸ＭＩＰ画像の対応する位置の画素に代入する。すなわち、ｚ軸ＭＩＰ画像は、ｚ軸上に並ぶ複数の画像から最大輝度の値を取り出して新たに作成した画像である。

なお、本実施例では、コンピュータの演算量を減少させるために３次元画像データを、２次元に投影してｚ軸ＭＩＰ画像を作成しているが、原理的には、３次元画像データのまま処理をすることももちろん好ましい。

５−２．ｚ軸ＭＩＰ画像を段階的にぼかす（スケールスペース変換）
次に、図１５のステップＳ１５−２に示すように、ｚ軸ＭＩＰ画像を段階的にぼかした画像を作成する。具体的には、ｚ軸ＭＩＰ画像に対して上述したスケールスペース変換を適用する。この際、σの値を段階的に変化させて段階的にぼかした複数の画像群を作成する。例えば、σ＝０．２μｍ、σ＝０．３μｍ・・・などのように段階的なスケールσを用いて複数の画像群を作成する。

例えば、σの値を、０．２μｍから５．０μｍまで、０．１μｍステップで変化させて、計４９段階の画像群を作成する。これらの複数段階の画像群を、「ステップぼかし画像群」と呼ぶ。

５−３．各スケールの画像について中心度の計算を行う
次に、図１５のステップＳ１５−３に示すように、上記ステップぼかし画像群に基づき、中心度の計算を行う。

まず、樹状突起を抽出するためには、境界線の存在度を中心度として用いる。中心度は、各画素毎に求める。そして、各画素毎にこの中心度の値を代入した画像を作成する。この画像を「中心度画像」と呼ぶ。

このような計算を、上記複数のステップぼかし画像群の各画像について行い、複数の「中心度画像」を作成する。

なお、本実施の形態において微分の演算がしばしば出てくるが、この微分は一般には、差分に置き換えて計算する。分点公式などを使うのが普通であり、本実施例でもそのような従来技術を利用することももちろん好ましい。

しかし、本実施例においては、畳み込み積分の形で書けるぼかし変換を行っているので、画像への微分演算が、畳み込みカーネルへの微分演算に帰着する。すなわち、最初に微分したカーネル関数を用意し、これを畳み込み積分することで、ぼかし変換と微分を同時に行っているのである。

この結果、本実施例によれば、コンピュータの演算量を減らすことができ、解析速度の向上を図ることができる。

５−４．中心度画像に対して、σ軸方向の最大中心度スケール投影を行う
次に、図１５のステップＳ１５−４に示すように、上記中心度画像に対して、σ軸方向での最大輝度の値を見つける。ここで、最大輝度とは、中心度の最大値を言う。そして、その最大輝度を与えるスケールσの値を取り出し、そのσの値を対応する画素に代入した「最大中心度スケール画像」を求める。

上述したように、σの値を変化させた複数のステップぼかし画像群（例えば２０枚）に対して、中心度画像がそれぞれ求められている。したがって、複数枚（例えば２０枚）の中心度画像群が得られるので、それらの画像郡から同じ位置の画素であって最大輝度（最大中心度）の画素を求め、その最大輝度の値を与えるスケールσを、画素に代入した新たな画像を形成する。

この画像は、いわばσ軸上で最大輝度（最大中心度）の画素の位置（すなわちσの値）を代入した画像である。この画像を「最大中心度スケール画像」と呼ぶ。

したがって、作成した投影画像（最大中心度スケール画像）の各点の値は、σ＝０．３μｍ、σ＝ ０．５μｍ、・・・のようになり、この値は、中心度が最も大きかったスケールσの位置を表す。

５−５．スネークを動かし、エネルギー極小となる位置を探す
次に、図１５のステップＳ１５−５に示すように、このようにして得られた最大中心度スケール画像中でスネークを動かし、エネルギー極小となる位置を探す。

さて、ここでは、スネークを複数のシリンダ図形を連結させて、太さが一様でない線図形としたものとして取り扱う。ここで、シリンダ図形とは、切頂円錐（ｔｒｕｎｃａｔｅｄ ｃｉｒｃｕｌａｒ ｃｏｎｅ）」である。すなわち、任意の円錐を平面で切断して、頂点を含む部分を除去した図形をシリンダ図形と呼ぶ。このように、スネークとは、シリンダ図形によって構成されたものとして扱うのが便利である。

上述した最大中心度スケール画像中に存在する輝度の谷線（ｃｌｅｆｔ ｌｉｎｅ）を検出する。これによりスネークの空間パラメータの内、ｘｙ成分のフィッティングと太さのフィッティングが完了する。

５−６．スネークの調整
次に、図１５のステップＳ１５−６に示すように、完成したスネークを元の三次元断層画像中に戻し、ｚ軸成分のフィッティングと太さの修正を行う。

５−７．スパイン頭部の検出
上述のようにスネークを構成し、適宜調整することによって、樹状突起の形態の検出ができた。次に、上記ステップ５−４で最大中心度スケール画像から、樹状突起の「樹状」の部分が得られるので、この樹状の部分を除いた残りの画像中にスパインがあるはずである。したがってこの残存部分の画像中から、図１５のステップＳ１５−７に示すようにスパイン頭部の検出を行う。

まず、上記５−４（図１５のステップＳ１５−４）の処理後、ヘステンソルの固有値の積を求めて、ピーク位置の抽出を行い、これをスパイン頭部（の中心位置）とする。

上述したように、このような処理は、「樹状」部分以外の領域に対して行う。

そして、この中心位置を含むような楕円体であって、スパイン頭部との誤差を最小にするような楕円体を構成（すなわち近似）し、この楕円体をスパイン頭部と認定する。

また、スパイン頭部の中心位置から樹状突起へ垂線を降ろし、スパインの柱部を構成する。この垂線（柱部）の長さを調べることによって、柱部の長さ（ｈ）を知ることができる。

このようにして得られたスパイン頭部（及び柱部）を先に抽出した樹状突起（すなわちスネーク）と合わせて神経細胞の形態の検出が完了する（ステップＳ１５−８）。

本実施例では、以上のようにして、神経細胞の形態の検出を行っているので、スパインの形状を含む神経細胞の形状を効率よく観察することができる。

５−８．変形例
なお、本実施例では、途中で二次元に投影しているが、これはコンピュータにおける処理時間を削減するためであり、理論的には、三次元のまま計算することももちろん好ましい。その場合は、例えば上記中心度画像は四次元（たて・横・高さ・スケール）の画像になる。
「実施例２」

６．実施例２：改良された神経細胞解析方法
６−１．背景
上述した実施例１では、樹状突起のトレース（４−１−２．）の際にスネーク（関数）を用いている（数２）が、計算量が過大になりがちであり、計算が長時間に渡ってしまう場合もある。特に、細かい構造が多数存在する場合に、それらを全て正確に求めるのには、膨大な時間がかかる場合もある。上記５−５．で述べたようにスネークを動かしてエネルギー最小の位置を探す作業は計算量が非常に大きなものとなる傾向にある。その結果、診断に支障が出る恐れも想定される。

そこで、より迅速に樹状突起のトレース等の演算を行うための改良された手法を以下実施例２として提案する。

６−２．原理
本実施例２で提案する方法は、スケールスペース法とfacet（切子）モデルとを組み合わせて、画像中から線状図形を検出する方法である。Facet モデルとは、ピクセルの大きさよりも微細な構造について論じるための枠組みであり、各ピクセルの内部を多項式で近似して、その集合として画像を捉えることを特徴とする。

二次元画像は、実数値関数

として表すことができる。画像f が十分に滑らかな場合は、定義域の任意の点の近傍で次式が成り立つ（Taylor の定理）。

ここで、標本化が十分な精度で行われている場合、すなわち、ピクセルの大きさが上記の展開が可能な近傍の大きさよりも十分に小さい場合は、各ピクセルの内部を次の多項式で近似することができる。

ただし、直感的な理解を助けるため、一階微分と二階微分を標準基底で表示し、それぞれgrad f、H と表した。ここで、さらに二階微分テンソル（ヘステンソル）の対角化基底を取って表示すると、上式（数７）はさらに簡略化されて、

と表すことができる。ここで、

は、ヘステンソルの二つの固有値のうち、大きいほうと小さいほうをそれぞれ表しており、

は対応する方向の微分係数を表す。今、注目しているピクセルの上を線状図形が通過しているとき、その図形の大きさを十分に捉えうるスケールのもとで、ヘステンソルの固有値は一方が大きな負の値を取り、他方が０に近い値を取る。この時、小さい方の固有値（負の固有値）に属する固有ベクトルは線状図形の法線方向を示し、法線方向の主軸と呼ばれる。これに対し、大きい方の固有値に属する固有ベクトルは接線方向を示すので、接線方向の主軸と呼ばれる。もし、ピクセルが図形の中心線上にあるならば、そこは尾根状になっているはずであるので、法線方向で極値を取らなければならない。すなわち、

の直線上で

を解いて、極値の位置

を得るが、この変位

が注目するピクセル内に存在しなくてはならない。したがって、定義域のすべての点について上記の変位を求め、それがピクセル内に含まれるか調べることで、その点が中心線上にありうるかどうか判別することができる。

６．３ 処理フロー
以下、具体的な処理の流れをフローチャートに基づき説明する。図１６には実施例２に係る動作を表すフローチャートが示されている。この処理フローは、樹状突起のトレースの処理フローであり、上述した４−１−２．の「樹状突起のトレース」で説明した処理に代わる処理であり、具体的には、５−１．から５−５．にいたる処理の変形例と考えることができる。

まず、ステップＳ１６−１においては、最適なスケールの選択を行う。スケールの選択方法については後に詳述する。

次に、ステップＳ１６−２において、ＭＩＰ（Max Intensity Projection）画像を作成する。

次に、ステップＳ１６−３においては、ＭＩＰ画像の各ピクセルについて、小さい方の固有値とそれに対応する固有ベクトルを求める。すなわち、

と法線方向の算出を行う。

ステップＳ１６−４においては、ＭＩＰ画像の各ピクセルについて、一階微分を求める。すなわち、

の算出を行う。

ステップＳ１６−５においては、ＭＩＰ画像の各ピクセルについて、中心度を求める。ここで、中心度とは、ピクセル内に「実施例１」で述べた極値が存在すれば

となり、存在しなければ０となる量である。

ステップＳ１６−６においては、中心度がゼロでない点で、その上を通る線状図形の直径を求める。直径は、以下のサブプロセスによって求めることができる。

（ステップａ）注目している点をＯとし、そこから法線方向の両側に等距離だけ進んで得られる二点をA、B とする。

（ステップｂ）点Ａ、Ｂでそれぞれベクトル

に対する境界度を求める。ここで、ベクトルｅに対する境界度βとは、その点における微分係数をgrad fとしたとき、β＝−ｅ・grad fで定義される量である。

（ステップｃ）求めた２つの境界度の積をとる。この積を

で表す。

（ステップｄ）点Ｏから点Ａ、Ｂまでの距離ｒを０．１?５．０μｍの範囲で動かし、それぞれの距離で境界度積

を求める。

（ステップｅ）境界度積

のピークのうち、中心Ｏに最も近いものの位置を直径の端点として確定する。この時、直径の値は、この端点から点Ｏまでの距離の二倍として得られる。

ステップＳ１６−７においては、上記ステップＳ１６−５で得られた中心度の画像を走査し、中心度が最大の点を選び出す。さらに、この点から１μｍ離れた点の中で最大の中心度を持つ点を選び出す。選び出したこれら２点を直線で結んで最初のスネークとする。これを初期スネーク（又は「初期のスネーク」）と呼ぶ。

ステップＳ１６−８においては、（初期）スネークの両端に、新しい端点を順次追加・接続し、スネークを成長させる。ただし、新しい端点は以下の条件により選出される。

（条件ａ）これまでの端点から新しい端点までの距離は１μｍ。

（条件ｂ）これまでの端点と新しい端点とを結んだ直線は、伸長前のそれに対して±３０度以内の角度をなす。

（条件ｃ）上記（条件ａ）、（条件ｂ）を満たす点の中で、最大の中心度を持つ点を新しい端点とする。新しい端点をスネークに加える。

（条件ｄ）端点をスネークに加えていき、スネーク（の端点）が画像の境界及び除外領域に到達したら、新しい端点の探索を停止する。画像中で新しい端点となる点を見つけられない場合に停止するのである。スネークの両方の端で探索が停止した場合に、新たな端点の探索は終了し、スネークの成長は完了し、最終的なスネークが作成される。

なお、このようにして見つかった端点は、最終的にはスネークのノード（節点）となり、一番最後に見つかった端点は、スネークの両端の端点となる。

ステップＳ１６−９においては、作成したスネークの各点（ノード、節点）で直径を取得し、スネーク内部の領域を確定する。この内部領域を画像の除外領域に追加する。

ステップＳ１６−１０においては、作成したスネークが３つ以上のノードからなり、かつ、全ノードの平均中心度が、これまでに確定した樹状突起のそれの１０％以上である場合は、この得られたスネークを樹状突起として確定する。

一方、この条件を満たさない場合は、ここで作成したスネークを破棄してステップＳ１６−７に処理が戻る。この操作は、画像からスネークが取り出せなくなるまで、もしくは最初に決めた回数だけ繰り返して終了する。

以上のようにして樹状突起のトレースを行う。

スケールの選択
上述したように、ステップＳ１６−１におけるスケールの選択方法についてその詳細を述べる。スケール範囲σ＝０．１μｍ?３．０μｍで、操作者が自由に値を選べるようにし、操作者が選択した値に対してステップＳ１６−２から、Ｓ１６−６までを実行する。

ステップＳ１６−６中の（ステップｅ）で得られた境界度積のピークの高さ（位置ではない）を各点に対してプロットして得られた画像を操作者に対して表示する。

操作者はこの画像を見、この得られた画像が適切に中心線を表していると判断した場合は、ステップＳ１６−７以下の処理が実行される。適切に表していないと判断した場合には、再び別のスケールの選択を行うため、ステップＳ１６−１の処理に移行する。

本実施例２によれば、スネークのフィッティングに代えて、スネークの端点を順次追加して成長させる処理を実行しているため、計算量を減らすことができ、迅速な診断に寄与するものである。

７．その他の変形例・応用例
（１）上述した実施の形態では、神経細胞の解析方法について説明した。しかしながら、この手法は、神経細胞のように線状、樹状の形態を有する他の細胞、組織、器官であればそのまま利用できることは言うまでもない。

例えば、毛細血管のように網状に形成されている組織に対しても、本実施の形態の手法を応用することが可能である。

特に、本実施の形態の手法によれば、線状の形態の上に、棘状の部位が発生するような場合も、その棘状部位を効率的に解析することが可能である。棘状部位とは典型的にはスパインが挙げられる。

（２）また、本実施の形態では、細胞の形態の３次元画像を入力として用いた。このような３次元画像は、共焦点顕微鏡を用いて容易に撮影することが可能であるが、他の手法・装置で３次元画像を得てもかまわない。例えば、２光子顕微鏡等その他の種々のイメージング装置を利用可能である。

（３）また、上記実施例では、コンピュータにプログラムをインストールし、神経細胞の３次元画像データを、このコンピュータで解析する例を説明したが、このプログラムは、コンピュータのハードディスク等に予めインストールしておくことが好ましい。また、プログラムは、ＣＤＲＯＭ等に記録しておき、所望のコンピュータに供給することが好適であるが、ネットワークを介してコンピュータにインストールすることも好ましい。

また、入力となる３次元画像データは、一般的な２次元の画像フォーマットの画像を複数枚用いる手法で実現し、コンピュータに供給することが好ましい。もちろん、特定の３次元データ専用のデータフォーマットを適宜作成することも好ましい。

この３次元画像データは、ＣＤＲＯＭや各種磁気ディスク、光ディスク、半導体記憶装置に一旦格納してから、コンピュータに供給することが好ましい。また、顕微鏡とコンピュータとを通信回線で接続して、３次元画像データを転送する構成を採用することも好ましい。

（４）また、上記実施例では、汎用的なコンピュータを用いて、神経細胞解析装置を実現したが、専用ハードウェアを用いて神経細胞解析装置を構築することも好ましい。

特に、顕微鏡と一体構造にする場合や、検査装置に組み込む場合には、汎用的なコンピュータにプログラムをインストールする構成よりも、専用の装置を構成した方が便利である場合もある。

Claims (25)

1. 線状の形態を有する細胞又は組織又は器官である解析対象物の３次元画像に基づき、前記対象物を解析する方法において、
   前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかしステップと、
   前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度ステップと、
   前記第１の中心度に基づき、前記対象物の線状の形態に近似した図形を求める線状形態フィッティングステップと、
   を含むことを特徴とする線状の形態を有する細胞又は組織又は器官を解析する方法。
2. 神経細胞の３次元画像に基づき、神経細胞の形態を解析する方法において、
   前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかしステップと、
   前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度ステップと、
   前記第１の中心度に基づき、神経細胞の樹状突起に近似した図形を求める樹状突起フィッティングステップと、
   を含むことを特徴とする神経細胞解析方法。
3. 請求項２記載の神経細胞解析方法において、
   前記第１のぼかしステップは、
   前記第１のスケールσの値を変化させてスケールスペース変換を行い、各スケールσに対する複数のぼかし画像群を生成することを特徴とする神経細胞解析方法。
4. 請求項３記載の神経細胞解析方法において、
   前記第１のぼかしステップは、
   前記第１のスケールσの値を０．２μｍ−５．０μｍの範囲で変化させてスケールスペース変換を行い、各スケールσに対する複数のぼかし画像群を生成することを特徴とする神経細胞解析方法。
5. 請求項２記載の神経細胞解析方法において、
   前記第１の中心度ステップは、境界面の存在度を中心度として求めることを特徴とする神経細胞解析方法。
6. 請求項２記載の神経細胞解析方法において、
   前記第１の中心度ステップは、はぎ取り法を用いて中心線を求めることを特徴とする神経細胞解析方法。
7. 請求項２記載の神経細胞解析方法において、
   前記樹状突起フィッティングステップは、前記求めた中心度に基づき、スネークを樹状突起にフィッティングさせることを特徴とする神経細胞解析方法。
8. 線状の形態を有する細胞又は組織又は器官である解析対象物の３次元画像に基づき、前記対象物を解析する装置において、
   前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかし手段と、
   前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度手段と、
   前記第１の中心度に基づき、前記対象物の線状の形態に近似した図形を求める線状形態フィッティング手段と、
   を含むことを特徴とする線状の形態を有する細胞又は組織又は器官を解析する装置。
9. コンピュータを、線状の形態を有する細胞又は組織又は器官である解析対象物の３次元画像に基づき、前記対象物を解析する装置として動作させるプログラムにおいて、前記コンピュータに、
   前記３次元画像に対して、所定の第１のスケールσに基づきスケールスペース変換を行う第１のぼかし手順と、
   前記スケールスペース変換後の前記３次元画像の各画素に対して、第１の中心度を求める第１の中心度手順と、
   前記第１の中心度に基づき、前記対象物の線状の形態に近似した図形を求める線状形態フィッティング手順と、
   を、実行させることを特徴とするプログラム。
10. 請求項１記載の線状の形態を有する細胞・組織又は器官を解析する方法において、
    前記３次元画像に対して、所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかしステップと、
    前記第２のぼかしステップによるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度ステップと、
    前記第２の中心度に基づき、前記線状の形態の上に存在する棘状部位の検出を行う棘状部位検出ステップと、
    を含むことを特徴とする線状の形態を有する細胞・組織又は器官を解析する方法。
11. 請求項２記載の神経細胞解析方法において、
    前記３次元画像に対して、第１のスケールσとは異なる所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかしステップと、
    前記第２のぼかしステップによるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度ステップと、
    前記第２の中心度に基づき、神経細胞のスパイン頭部の中心位置を検出するスパイン頭部検出ステップと、
    を含むことを特徴とする神経細胞解析方法。
12. 請求項１１記載の神経細胞解析方法において、
    前記第２のスケールσの値を０．１μｍ−０．３μｍの範囲で変化させてスケールスペース変換を行い、各スケールσに対する複数のぼかし画像群を生成することを特徴とする神経細胞解析方法。
13. 請求項１１記載の神経細胞解析方法において、
    前記スパイン頭部検出ステップは、前記第２の中心度の絶対値がピークを取る位置を、スパインの頭部の中心位置と判断することを特徴とする神経細胞解析方法。
14. 請求項１１記載の神経細胞解析方法において、
    前記スパイン頭部検出ステップは、前記樹状突起から所定の距離以内の位置を、スパインの頭部の中心位置と判断することを特徴とする神経細胞解析方法。
15. 請求項１１記載の神経細胞解析方法において、
    前記求めたスパイン頭部の中心位置から、最も近い前記樹状突起へ、垂線を降ろし、その垂線をスパインの柱部として求めるスパイン柱部検出ステップ、
    を含むことを特徴とする神経細胞解析方法。
16. 請求項１１記載の神経細胞解析方法において、
    前記スパイン検出ステップは、ヘステンソルの固有値が全て負である領域をスパインの頭部の存在領域として求めることを特徴とする神経細胞解析方法。
17. 請求項１６記載の神経細胞解析方法において、
    前記スパイン検出ステップは、前記求めた存在領域を、楕円体で近似し、スパイン頭部の近似図形を求めることを特徴とする神経細胞解析方法。
18. 請求項９記載の線状の形態を有する細胞・組織又は器官を解析する装置において、
    前記３次元画像に対して、所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかし手段と、
    前記第２のぼかし手段によるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度手段と、
    前記第２の中心度に基づき、前記線状の形態の上に存在する棘状部位の検出を行う棘状部位検出手段と、
    を含むことを特徴とする線状の形態を有する細胞・組織又は器官を解析する装置。
19. 請求項１０記載の、コンピュータを、線状の形態を有する細胞・組織又は器官を解析する装置として動作させるプログラムにおいて、
    前記３次元画像に対して、所定の第２のスケールσに基づきスケールスペース変換を行う第２のぼかし手順と、
    前記第２のぼかし手順によるスケールスペース変換後の前記３次元画像の各画素に対して、ヘステンソルの固有値の積を第２の中心度として求める第２の中心度手順と、
    前記第２の中心度に基づき、前記線状の形態の上に存在する棘状部位の検出を行う棘状部位検出手順と、
    を実行させることを特徴とするプログラム。
20. 請求項２記載の神経細胞解析方法において、
    前記樹状突起フィッティングステップは、
    前記求めた中心度に基づき、中心度が最大の点を選択する第１選択ステップと、
    前記中心度が最大の点から所定距離離れた点の中で、前記中心度が最大の点を選択する第２選択ステップと、
    前記選択した２個の点を結び、初期スネークを構成するステップと、
    前記構成した初期スネークの端点から所定距離離れた点の中から、最大の中心度を持つ点を選択し、既存の端点と結び、前記スネークを成長させる成長ステップと、
    を含むことを特徴とする神経細胞解析方法。
21. 請求項２０記載の神経細胞解析方法において、
    前記成長ステップにおいて、
    選択する点は、前記既存の端点との距離が所定の距離にある点であることを特徴とする神経細胞解析方法。
22. 請求項２０又は２１記載の神経細胞解析方法において、
    前記成長ステップにおける選択する点を見つけられる場合、前記成長ステップにおける成長を続行し、前記選択する点を見つけられない場合は、前記成長ステップにおける成長を完了する制御ステップ、
    を含むことを特徴とする神経細胞解析方法。
23. 請求項２２記載の神経細胞解析方法において、
    前記成長ステップにおける選択する点を見つけられる場合とは、前記選択する点を、画像中において既にスネークの内部と認定されている領域を除いた領域から見つけられる場合であることを特徴とする神経細胞解析方法。
24. 請求項２０、２１、２２、又は２３記載の神経細胞解析方法において、
    前記成長させたスネークの各ノードで直径を取得し、スネーク内部の領域を確定させる確定ステップ、
    を含むことを特徴とする神経細胞解析方法。
25. 請求項２０、２１、２２、２３、又は２４記載の神経細胞解析方法において、
    前記成長させたスネークが、
    （１）３個以上のノードを有すること、
    （２）全ノードの平均中心度が、所定のそれまでに確定した樹状突起の平均中心度の所定の比率以上であること、
    の２条件をともに満たす場合に樹状突起として確定し、いずれかの条件を満たさない場合には、そのスネークを破棄し、前記第１選択ステップに戻り、中心度が最大の点を選択する動作を再び開始する繰り返し制御ステップ
    を含むことを特徴とする神経細胞解析方法。

Images