JP5935146B2 - 眼科疾患の画像解析装置、眼科用画像解析方法及び眼科用画像解析プログラム - Google Patents

Description

本発明は、眼科疾患の画像解析技術に関し、特に、強度近視の画像解析技術に関する。

眼科疾患でポピュラーな近視の中に、眼球後極の組織が変性して眼軸が伸び、強度な近視の症状を呈し、場合により失明に至る重篤な疾患「強度近視」が知られている。強度近視は、網膜の変性や萎縮などが原因で起こり、眼球が楕円形に伸び、眼軸長が正視より３．５ミリ以上長い。全国で４０歳以上の人口約６７００万人のうち、３６０万人ほどが強度近視とみられている。網膜に負担がかかり、飛蚊症、網膜剥離、緑内障などの病気をひきおこし、視野が狭くなり、視力低下や失明にいたる可能性がある。

現状では本疾患に対する根治術は特許文献１など研究段階で、薬物療法や外科治療により病態の症状や進行を抑える対症療法しかなく、早期発見が重要となる。しかし、本疾患は２０年以上の長期に渡るタイムスパンで病態が徐々に進行し、早期の段階では眼軸長の延長や視力の低下を呈さないことが多いため、眼底カメラなど現状の眼科診断機器では発見が遅れがちになるという問題があった。

近視は、屈折矯正するために必要なレンズの度数Ｄ（ディオプターあるいはディオプトリ、Ｄ＝１（ｍ）/焦点距離（ｍ））の値を基準に判定される。正視の場合は、Ｄ値は０となり、遠視では正値となる。Ｄが負値の場合が近視で、−６Ｄ以下が強度近視である。ただし、強度近視の初期段階では、Ｄ値が異常値を示さないことがあるため、眼球が前後方向に伸びる特徴を基準にする方が早期に診断できる。

本疾患に対する従来の診断手法として、上記の屈折検査装置（特許文献２「眼屈折計」、キャノンなど）による視力測定（特許文献３「眼鏡、コンタクトレンズ度数決定システム」、ビジョンメガネ）のほかに、超音波による眼軸測定（特許文献４「眼科用超音波診断装置」、キャノン）、光干渉による眼軸測定（特許文献５「生体眼の寸法測定装置」、トプコン）、眼底カメラによる網膜合併症解析（特許文献６「眼底疾患の解析装置」、トプコン）、また緑内障との合併症の診断に液体パルスによる眼圧測定（特許文献７「非接触式眼圧計」、トプコン）などがある。

特表２００７−５３６２７２号公報 特公平５−０５４７７４号公報 特許第４０１４４３８号公報 特公平３−０５８７２６号公報 特許第２７２３９６７号公報 特許第３４７９７８８号公報 特公平４−０３０２９４号公報 特開２００７−１２１２５９号公報

Brain J.Curtin, "The Posterior Staphyloma of Pathlogic Myopia", Transaction of American Opthalmology Society, Vol. LXXXV, pp.67-86,(1997).

非特許文献１に記載に強度近視の病態については、１９７７年にＢ.Ｊ.Ｃｕｒｔｉｎが発表した眼底カメラ像に基づく１０タイプ分類が知られている。

また、特許文献２から７までに示される従来のその他の眼科診断機器を用いて、眼軸延長や視力低下といった「強度近視」の典型的症状の診断は可能であるが、眼球後極の変性・肥大といった病変については、眼球前方からアプローチする従来の検査手段では計測困難であった。

また、画像診断を導入して変性・肥大を診断する場合、いずれかの方向に非対称性があるという定性的判定は読影医間でバラツキを生じさせやすいという問題もある。

本発明の第１の目的は、眼軸延長や視力低下といった「強度近視」の典型的症状を呈する前段階で、本疾患の早期兆候である眼球後極の変性・肥大（後部ぶどう腫とよばれる）の度合いを非侵襲かつ定量的に解析することである。そして、早期に「強度近視」の診断を確定し、病態の進行を抑止できるようにすることである。

また、本発明の第２の目的は、画像診断にあたっては読影医の主観に依存する判定ができるだけ生じないように解析の自動化と定量化を実現することである。

本発明は、頭部疾患の診断ではポピュラーな３次元トモグラフィー技術を眼科分野に導入し、主として頭部のＭＲＩ(核磁気共鳴画像、Magnetic Resonance Imaging)により眼球部の３次元画像を取得し、重症化しやすいかといった本疾患の予後と相関のある眼球後極の形態特徴パラメータを非侵襲的、定量的かつ自動的に解析する手法を提案するものである。

例えば、頭部のＭＲＩ画像より特許3049021号（「画像処理装置」、五大）等に記載の技術を用いて切り出した左右の眼球の３次元ＭＲボクセル画像を、特許44979651号（「医用画像表示装置」、日立メディコ）等の技術を用いて上下左右前後の６方向でボリュームレンダリングした６フレームの投影画像を準備し、眼球の解析を行う。

続いて、上下左右方向の４枚の各投影画像に対して眼軸を設定し、重心から後極部に向かう放射状のエリアで左右方向（耳鼻方向）または上下方向の幾何学的対称性と、後部曲面の尖鋭度を算出する。更に、オプションとして後部投影画像より後部の突起（ぶどう腫）の個数が単一か複数かといった解析を行う。

眼軸の設定は自動抽出を基本とし、場合により対話的手段により補正を加える。これ以外の算出処理は全て自動化し、操作者間で算出される各パラメータにバラツキが生じないようにする。これらの眼軸長、左右対称性、上下対称性、尖鋭度、突起数の４＋オプション１項目を基に、強度近視の病態を分類し、予後判定や治療方針の策定に活用する。

本発明の一観点によれば、眼球を含む頭部の３次元ボクセル画像を取得する画像入力部と、入力した前記３次元ボクセル画像より、左右眼球別に眼球領域の少なくともいずれか一方の眼の３次元眼球ボクセル画像を抽出するボクセル画像抽出部と、抽出した３次元眼球ボクセル画像を、少なくとも上下いずれかの方向に投影した第１の投影画像を作成する投影画像算出部と、前記第１の投影画像に対して眼球前極部頂点と眼球後極部頂点とを結ぶ眼軸を設定し、少なくとも眼球前極部頂点と眼球後極部頂点との間の眼軸長と、眼軸を境界に後極側を左右に分割した領域をもとに左右方向の非対称性の度合いと、を画像解析のパラメータとして算出する投影画像解析部と、を備えることを特徴とする眼科疾患の画像解析装置が提供される。

この発明によれば、眼軸長の長さ又は左右方向の非対称性の度合いに基づいて、強度近視を判定することができる。パラメータは強度近視を特定するためのパラメータを含む。

前記３次元ボクセル画像は、ＭＲＩを用いて取得したものであることを特徴とする。
前記投影画像算出部は、前記３次元眼球ボクセル画像を、少なくとも左右いずれかの方向に投影した第２の投影画像を作成し、前記画像解析部は、前記第２の投影画像に対して眼軸を設定し、少なくとも眼球前極部頂点と眼球後極部頂点との間の眼軸長と眼軸を境界に後極側を上下に分割した領域をもとに上下方向の非対称性の度合いを、画像解析のパラメータとして算出することを特徴とする。これにより、上下方向の非対称性を算出することができる。

また、前記画像解析部は、前記第１又は第２の投影画像のうちの少なくともいずれかに対して、後極の円弧部曲率の尖鋭度合いを、画像解析のパラメータとして算出することを特徴とする。これにより、後極尖鋭度の算出をすることができる。

また、前記画像解析部は、前記第１又は第２の少なくともいずれか一方の投影画像について、左右又は上下の少なくともいずれか一方の組の前記パラメータを平均した値を算出結果とすることを特徴とする。これにより、非対称性、後極尖鋭度の平均値を算出することができる。

また、前記投影画像算出部は、前記３次元眼球ボクセル画像に対して、後方向に投影した第３の投影画像を作成し、前記画像解析部は、前記第３の投影画像に対して、後極突出部に対応するピーク点の個数を、画像解析のパラメータとして算出することを特徴とする。これにより、後極突出数の算出をすることができる。

また、前記画像解析部は、前記第３の投影画像に対して、眼球部の輪郭線抽出を行い、更に眼球部の重心点を設定し、前記算出された各ピーク点の位置を、画像解析のパラメータとして、前記重心点を基準とした極座標である重心からの距離と重心からピーク点方向の角度で算出することを特徴とする。これにより、後極突出数の位置を算出することができる。

また、前記画像解析部は、前記投影画像に対して眼球部の輪郭線抽出を行い、前記眼球前極部頂点および前記眼球後極部頂点は前記眼球部の輪郭線上に位置し、かつ前記眼球前極部頂点と前記眼球後極部頂点とを結ぶ直線が前記直線の幾つかの頂点より鉛直方向に左右または上下方向に延長した前記輪郭線との交点までの距離が均等になるように、前記眼球前極部頂点および前記眼球後極部頂点を設定することにより、前記眼軸を自動設定することを特徴とする。これにより、眼軸の自動設定をすることができる。

或いは、前記画像解析部は、前記投影画像に対して眼球部の輪郭線抽出を行い、更に眼球部の重心点を設定し、前記眼球前極部頂点および前記眼球後極部頂点は前記眼球部の輪郭線上に位置し、かつ前記眼球前極部頂点と前記重心点と前記眼球後極部頂点との３点を結ぶ折れ線ができるだけ直線になるように、前記眼球前極部頂点、前記眼球後極部頂点および前記重心点を設定することにより、前記眼軸を自動設定するようにしても良い。

また、前記画像解析部は、前記第１又は第２の投影画像に対して眼球部の輪郭線抽出を行い、前記眼球前極部頂点または前記眼球後極部頂点のいずれかに対して、画面との対話形式により設定する対話型頂点入力部により設定された前記眼球前極部頂点または前記眼球後極部頂点に対して、前記眼球前極部頂点と前記眼球後極部頂点とを結ぶ直線の延長線上で前記輪郭線との交点に自動的に補正する頂点補正部を備えていることを特徴とする。これにより、眼軸の手動設定をすることができる。

また、前記画像解析部は、前記第１又は第２の投影画像に対して眼球部の輪郭線抽出を行い、更に眼球部に計測用重心点を設定し、計測用重心点より前記眼球後極部頂点を結ぶ基準線に対して、計測用重心点を中心に＋第１の角度および−第１の角度の方向の直線Ａと直線Ｂを設定し、前記直線Ａと前記基準線に囲まれた前記輪郭線の内部の面積Ａと、前記直線Ｂと前記基準線に囲まれた前記輪郭線の内部の面積Ｂを算出し、面積Ａと面積Ｂとの比率を前記非対称性の度合いとして与えることを特徴とする。これにより、非対称性の算出を行うことができる。
尚、第１の角度は、例えば、４０度以上９０度未満である。

前記画像解析部は、前記第１又は第２の投影画像に対して眼球部の輪郭線抽出を行い、更に眼球部に計測用重心点を設定し、計測用重心点より前記眼球後極部頂点を結ぶ基準線に対して、計測用重心点を中心に＋第２の角度および−第２の角度の方向の直線Ａと直線Ｂを設定し、前記直線Ａと前記輪郭線との交点Ａと前記眼球後極部頂点と前記直線Ｂと前記輪郭線との交点Ｂとの３頂点のなす角を前記尖鋭度合いとして与えることを特徴とする。これにより、後極尖鋭度の算出をすることができる。尚、第２の角度は、例えば、１０度から４０度程度である。

或いは、前記画像解析部は、計測用重心点より前記眼球後極部頂点を結ぶ基準線に対して、計測用重心点を中心に＋第２の角度および−第２の角度の方向の直線Ｃと直線Ｄを追加設定し、前記直線Ｃと前記基準線に囲まれた前記輪郭線の内部の面積Ｃと、前記直線Ｄと前記基準線に囲まれた前記輪郭線の内部の面積Ｄを追加算出し、面積Ｃと面積Ｄとの比率と前記面積Ａと面積Ｂとの比率のいずれか大きい値を前記非対称性の度合いとするようにしても良い。

また、前記画像解析部は、前記計測用重心点として前記眼球前極部頂点と前記眼球後極部頂点とを結ぶ眼軸上に、前記眼球後極部頂点より所定の距離だけ離した点を与えることを特徴とする。これにより、計測用重心点の設定をすることができる。

例えば、前記画像解析部は、前記所定の距離として正常者の眼球の平均的な眼軸長の１／２（例えば、８７画素、１３．４ｍｍ）を設定するようにしても良い。

本発明の他の観点によれば、眼球を含む頭部の３次元ボクセル画像を取得する画像入力ステップと、入力した前記３次元ボクセル画像より、左右眼球別に眼球領域の少なくともいずれか一方の眼の３次元眼球ボクセル画像を抽出するボクセル画像抽出ステップと、抽出した３次元眼球ボクセル画像を、少なくとも上下いずれかの方向に投影した第１の投影画像を作成する投影画像算出ステップと、前記第１の投影画像に対して眼球前極部頂点と眼球後極部頂点とを結ぶ眼軸を設定し、少なくとも眼球前極部頂点と眼球後極部頂点との間の眼軸長と、眼軸を境界に後極側を左右に分割した領域をもとに左右方向の非対称性の度合いと、を画像解析のパラメータとして算出する投影画像解析ステップと、を備えることを特徴とする眼科用画像解析方法が提供される。

又、本発明は、コンピュータに、上記に記載の眼科用画像解析方法を実行させるための眼科用画像解析プログラムであっても良く、当該プログラムを記録するコンピュータ読み取り可能な記憶媒体であっても良い。

本発明では、ＭＲＩを用いて強度近視の３次元的な眼球特徴を非侵襲的に計測することにより、臨床的な予後と相関が高そうな形状特徴パラメータに基づく、新規な病態分類手法を提案する。６方向からレンダリングした３次元ＭＲ投影画像を基に、合併症の併発と相関が高そうな、眼軸長、眼球後部の左右および上下の対称性、眼球後部の尖鋭度を自動的に算出することができ、これらを基に病態分類できる。

本発明によれば、患者に負担をかけずに「強度近視」の早期診断に寄与でき、読影医の主観に依存しない客観性のある診断データを提供することができる。

また、安価な汎用パーソナルコンピュータのソフトウェアで画像処理を実現することができ、ＧＰＵ（Graphics Processing Unit）など特殊なグラフィックスボードの増設も必要とないため、安価な診断サービスを提供できる。

強度近視診断の大まかな診断処理の流れを示すフローチャート図である。 本実施の形態による、頭部ＭＲＩ撮影データを用いた眼科疾患の画像解析技術を利用した強度近視診断手法の概要を示すフローチャート図である。 図２Ｂ（ａ）は、上部から撮影したＭＲＩ投影画像の例を示し、図２Ｂ（ｂ）は左から見た眼球の投影画像、図２Ｂ（ｃ）は上方から見た眼球の投影画像、図２Ｂ（ｄ）は後方から見た眼球の投影画像である。 図３（ａ）は、本発明の一実施の形態による画像解析装置の一構成例を示す機能ブロック図であり、図３（ｂ）は、投影画像解析部の構成例を示す機能ブロック図である。 本実施の形態によるＭＲＩにより取得した画像を用いた強度近視の計測派多メータを示す図であり、例えば図２Ｂの（ｂ）から（ｄ）までの画像より求めるが、詳細な求め方は一例である。 本実施の形態による、ＭＲＩによる強度近視診断手法の一例を示すフローチャート図である。 前後左右上下の各画像（３００×３００画素）であり、前後の投影画像は解析に必須のものではないが、参照用に取得している。 ＭＲＩによる強度近視の計測パラメータの詳細を示す図である。 ＭＲＩによる強度近視の計測パラメータの詳細を示す図である。 ステップＳ３１−１の、上下方向又は左右方向の投影画像の解析処理を示す図である。 ステップＳ３１−１の、上下方向又は左右方向の投影画像の解析処理を示す図である。 ステップＳ３１−１の、上下方向又は左右方向の投影画像の解析処理を示す図である。 図５のステップＳ５１からの、後方向投影画像処理の様子を示す図である。 図１２（ａ）は、図１１（ａ）に続く処理を示す図である。図１２（ｂ）〜（ｅ）は、処理の様子を示す図である。 画像の二値化処理と輪郭抽出処理（Ｓ３２）の詳細について説明する図である。 図１４（ａ）は、近傍領域の統合処理の流れを示すフローチャート図であり、図１４（ｂ）〜（ｄ）は、処理の様子を示す図である。 眼軸の自動設定のアルゴリズムの概要を示す図である。 眼軸の自動設定処理（ステップＳ３３）の詳細を示す図である。 図１５に続く処理の流れを示す図である。 後極側エッジ点の補正と眼軸の決定処理とその様子を示す図である。 眼軸の修正処理の流れを示す図と、その様子を示す図である。 後極輪郭部におけるエッジ点の作成処理（Ｓ３８）の詳細とその様子を示す図である。 図２０（ａ）は、平滑化フィルタ処理（図１１（ａ）：Ｓ５２）の流れを示すフローチャート図である。図２０（ｂ）は、平滑化フィルタの例を示す図である。 投影画像に対する濃淡ピーク点の抽出処理（図１１（ｂ）：Ｓ５３参照）の流れを示すフローチャート図である。 図２１の処理で抽出できたピーク点の周辺画素のクリア処理の流れを示すフローチャート図である（図１１のステップＳ５４）。 原画像に対する重心点の算出処理を示すフローチャート図である（図１２（ａ）：ステップＳ５５）。 画素の二値化処理（ステップＳ３２）の詳細を示すフローチャート図である。 図２５（ａ）は、近傍画素の連結処理、領域抽出処理の流れを示すフローチャート図である（ステップＳ３２−３、図１３（ｄ）参照）。図２５（ｂ）は、８近傍画素のうちの４つの画素を示す図である。 図２６（ａ）は、近傍領域の統合処理の流れを示すフローチャート図である（図１４、ステップＳ３２−４、図１４（ｂ）参照）。図２６（ｂ）は、８近傍画素の例を示す図である。 図２６に続く処理の流れを示す図である。 周辺のごみのような微小面積を削除して、眼軸の決定における間違いを低減するための、微小面積領域の削除処理（ステップＳ３２−５）を示すフローチャート図である（図１４（ｃ）参照）。 輪郭抽出処理の流れを示すフローチャート図である（図１４（ｄ）参照）。 前極側および後極側エッジ点と重心点の算出処理の流れを示すフローチャート図である（Ｓ３３−１：図１５Ｂ）。 前極側エッジ点の補正と決定処理の流れを示すフローチャート図である（Ｓ３３−３、図１６参照）。 図３１のステップＳ３３−３−１の処理である前後エッジ点と重心点方向の距離差算出処理の流れを示すフローチャート図である（図１６、ステップＳ３３−４）。 重心点の補正と決定処理（Ｓ３３−４）の流れを示すフローチャート図である（図１６（ｃ）参照）。 後極側エッジ点の補正と眼軸の決定処理の流れを示すフローチャート図である（図１６の眼軸の自動設定処理Ｓ３３）。 ３頂点の直線性の判定処理の流れを示すフローチャート図である（図１８参照）。 手動設定エッジ点の輪郭線上への自動補正処理（Ｓ３５−３）を示すフローチャート図である。 既に決定された眼軸と直交する方向に９０度および−９０度エッジ点の追加を行う処理（ステップＳ３８−１、図１９（ａ）参照）の流れを示すフローチャート図である。 θ_ｗおよび−θ_ｗのエッジ点（±第１の角度）の追加処理の流れを示すフローチャート図である（ステップＳ３８−２、図１９（ａ）、（ｃ）参照）。 θ_ｎおよび−θ_ｎ度のエッジ点（±第２の角度）、の追加処理の流れを示すフローチャート図である（ステップＳ３８−３、図１９（ａ）、（ｄ）参照）。 （ａ）は左方向投影画像の場合の円弧面積の算出方法を示す図であり、（ｂ）は（ａ）で求めた面積の比率を示す図である。 図４０で説明した後極部における円弧面積の算出処理例を示すフローチャート図である。 （ａ）は、後極部における円弧面積の算出方法（狭領域）を示す図であり、（ｂ）は面積比を示す図である。 図４２で説明した後極部における円弧面積の算出処理例を示すフローチャート図である。 後極部３エッジ点（先鋭度）の成す角ηの算出方法を説明する図である。 解析対象の６方向投影画像例であり、左右、上下、前後の６方向からの投影画像が得られていることがわかる。 解析における解析パラメータの設定画面の例を示す図である。 解析結果の一例を示す図であり、解析処理に基づく解析結果を一画面で見る図である。 図４５に対応する図であり、解析時の補助線設定状況の確認画面の一例を示す図である。 強度近視の場合の解析対象の６方向投影画像の例を示す図である。 解析パラメータの設定画面例を示す図である。 図５０と同じ設定において、解析を行った結果の例を示す図である。 解析時の補助線設定状況の確認画面を示す図である。 解析対象の６方向投影画像（軽症）の６方向投影画像の例を示す図である。 左眼を解析対称とした解析パラメータの設定画面である。 図５４に基づく解析結果を示す図である。 解析時の補助線設定状況の確認画面の例を示す図である。 解析時の補助線設定状況の確認画面であって、眼軸の修正（この例では手動修正）を付加した様子を示す図である。 解析パタメータの設定画面であり、例えば図５４に示す設定画面から設定値が変更されていないが、設定値を画面で確認することができる。 修正後の再解析結果の一例を示す図である。

本発明の実施の形態における強度近視の検査技術は、３次元トモグラフィー技術を用いて、眼球を含む複数の視点からの投影画像を取得し、この投影画像に基づいて、強度近視の診断を行うことを特徴とするものである。

以下に、本発明の実施の形態による頭部ＭＲＩ撮影データを用いた眼科疾患の画像解析技術について、強度近視を例にして詳細に説明する。

まず始めに、強度近視診断の大まかな診断処理の流れについて、図１のフローチャート図を参照しながら説明を行う。まず、ステップＳ１において、一般的な眼科診断機器による診断を行う。一般的な眼科診断とは、視力検査、屈折率検査、網膜の病変検査、視野の検査などである。ステップＳ２において、診断は確定したか否かを判定する。診断が確定していれば、治療方針も確定する（ステップＳ８へ）。

ステップＳ２でＮｏの場合には、頭部ＭＲＩ撮影及び眼球部位のレンダリング（左右眼球別々に）を行う。ステップＳ４において、上下左右前後の６方向の頭部の投影画像データを左右眼球別に生成する。次いで、２次元（２Ｄ）投影画像に対する画像解析処理を行う。ステップＳ６において、形状パラメータおよび強度近視タイプの分類処理を行う。ステップＳ６の結果より、治療方針の検討を行い（ステップＳ７）、ステップＳ８に進み治療方針が確定する。

図２Ａは、本実施の形態による、頭部ＭＲＩ撮影データを用いた眼科疾患の画像解析技術を利用した強度近視診断手法の概要を示すフローチャート図である。図２Ｂ（ａ）は、人体の頭部を撮影したＭＲＩより眼球部位をスライスした断層画像の例を示し、更に頭部ＭＲＩ撮影データより片側の眼球部を抽出した３次元ボクセル画像に対して、図２Ｂ（ｂ）は左から見た眼球の投影画像、図２Ｂ（ｃ）は上方から見た眼球の投影画像、図２Ｂ（ｄ）は後方から見た眼球の投影画像である。

図２Ａに示すように、まず、ステップＳ１１で、頭部ＭＲＩ撮影を行い、眼球部位のボクセルデータを左右眼球別に抽出する（特許第３０４９０２１号公報に記載の方法を用いることができる）。

次いで、ステップＳ２１において、左右別に、各眼球ボクセルデータに対し、上下方向または左右方向の２次元投影像をボリュームレンダリング（上下方向または左右方向に各１枚または２枚）して、投影画像を作成する（特許第４４９７９６５号公報に記載の方法を用いることができる）。

次いで、ステップＳ３１において、左右眼球別に２〜４枚の上下方向または左右方向の各投影画像に対し眼軸を設定する。任意に、後極の形状特徴を投影画像により計測する。

さらに、ステップＳ５１おいて、左右眼球別に各眼球ボクセルデータに対し、後方向の２次元投影像をレンダリングして投影画像（各１枚）を作成する。この投影画像を用いて、任意に、突出数を計測する。

図２Ｂ（ａ）は、健常者の頭部ＭＲＩによる断層画像である。図２Ｂ（ｂ）の左から見た眼球の投影画像、図２Ｂ（ｃ）の上方から見た眼球の投影画像、図２Ｂ（ｄ）の後方から見た眼球の投影画像を用いている。

図３（ａ）は、本発明の一実施の形態による画像解析装置の一構成例を示す機能ブロック図であり、図３（ｂ）は、投影画像解析部の構成例を示す機能ブロック図である。

図３（ａ）に示す画像解析装置は、ＭＲＩにより取得した画像を入力する画像入力部１と、画像からボクセルデータを抽出するボクセル画像抽出部３と、ボクセルデータから投影画像を算出する投影画像算出部５と、この投影画像の画像解析を行う投影画像解析部７と、解析結果を出力する出力部１１と、を有している。図３（ｂ）に示すように、画像解析部７は、非対称性を算出する非対称性算出部７−１と、尖鋭度合いを算出する尖鋭度合い算出部７−２と、後極突出部のピークを算出する後極突出部ピーク算出部７−３と、眼軸を自動で設定する眼軸自動設定部７−４と、頂点を補正する頂点補正部７−５と、を有している。

図４は、本実施の形態によるＭＲＩにより取得した投影画像を用いた強度近視の計測パラメータを示す図であり、例えば図２Ｂの（ｂ）から（ｄ）までの投影画像より求める。図４（ａ）は、眼軸に沿った眼の長さである眼軸長Ｌ１を上下・左右の投影画像より求める方法を示す図であり、例えば、眼軸長Ｌ１を、上下・左右方向の４つの投影画像による計測結果の平均値で判定する。図２Ｂ（ｂ）の投影画像でいえば、白で表示されている眼球部分の中で左側の小さな半円領域（水晶体）の左端部分から、眼球の右端部分まで水平方向（図２Ｂ（ｂ）の場合は若干右上向きに）に伸ばした線の長さを求める。図４（ｂ）は、後極先鋭度θを求める方法を示す図であり、Ｃ１を眼軸として計測用重心より後極方向に±θｎの角度で伸ばした後極輪郭線との２つの交点と眼球後極部頂点とを結んだ線Ｌ_２、Ｌ_３のなす角度θ（後極尖鋭度θ）を、上下・左右方向の４つの投影画像による計測結果の平均値で判定する。図４（ｃ）は、左右対称性を求める方法を示す図であり、例えば、左右対称性を、上下方向の２つの投影画像による計測結果の平均値で判定する。具体的には、図４（ｃ）では、眼軸Ｃ１により区画された耳側領域（右目の場合：０〜＋θ_Ｗ、左目の場合：０〜−θ_Ｗ）の後極の面積Ｅと鼻側領域（右目の場合：０〜−θ_Ｗ、左目の場合：０〜＋θ_Ｗ）の後極の面積Ｎとの比を求める。実際には、後述するように、θ_Ｗより小さいθ_ｎを用いて面積Ｅと面積Ｎとの比を同様に求め、いずれか面積の比率の非対称度が大きい方を採用する。尚、この左右対称性は重症度との相関が高く重要なパラメータである。図４（ｄ）は、上下対称性を求める方法を示す図であり、例えば、上下対称性を眼軸Ｃ１により区画された上側領域Ｔ（０〜＋θ_Ｗ）と下側領域Ｂ（０〜−θ_Ｗ）との比Ｔ／Ｂ面積比率を求め、上下の値の平均値を求める。実際には、後述するように、θ_Ｗより小さいθ_ｎを用いて面積Ｔと面積Ｂとの比を同様に求め、いずれか面積の比率が大きい方を採用する。図４（ｅ）は、後極突出数を後方向の投影画像より求めるものであり、後部ぶどう腫の数が複数個あるか否かを求める。

表１は、ＭＲＩによる強度近視のタイプ分類の一例を示すものである。表１に示すように、４桁タイプの分類コード（００００〜２２２１までの全５４通り）は、１の位が眼軸長であり、正常が“０”、異常が“１”である。１０の位は、左右方向の対称性（上下の投影画像で後極左右方向θ_ｗ度円弧の面積比で判定する）であり、０：対称、１：鼻寄り（面積比９０％未満）、２：耳寄り（面積比１１０％以上）である。１００の位は、上下方向の対称性であり、左右像で後極上下方向θ_ｗ度円弧の面積比で判定する。０：対称、１：下寄り（９０％未満）、２：上寄り（１１０％以上）とする。１０００の位は、後極尖鋭タイプであり、上下左右の投影画像で後極先端中心の成す角で判定する。０：中間、１：紡錘型（１４０度未満）、２：樽型（１５０度以上）である。この４桁のコードにより、強度近視に関する病態等を分類して整理することができる。

図５は、本実施の形態による、ＭＲＩによる強度近視診断手法の一例を示す図である。図６は、眼球の前後左右上下方向の各投影画像（３００×３００画素）であり、前後方向の投影画像は解析に必須のものではないが、参照用に取得している。

まず、ステップＳ１１において、頭部ＭＲＩ撮影を行う（眼球部位のボクセルデータを左右眼球別に抽出する）。その処理アルゴリズムとしては、例えば、特許第３０４９０２１号に記載の方法を用いることができる。左右眼球別に各眼球ボクセルデータに対し、上下、左右、前後６方向の２次元投影像をレンダリングする（全６枚）。その処理アルゴリズムとしては、例えば、特許第４４９７９６５号に記載の方法を用いることができる。

より詳細には、上方向・下方向・左方向・右方向の４投影画像に対しては、後述する手法により眼軸を自動または手動で設定し、パラメータを算出する。後方向の投影画像については、後部ぶどう腫の突出の個数を自動または目視で判定するためのもので、複数個存在する場合に、後述する左右対称性または上下対称性の結果が適切に判定されない場合があるため、その結果を補足するために使用する。後部ぶどう腫の突出の個数を自動抽出する手法としては、上記特許文献８に記載の核医学画像よりピーク点を自動抽出するアルゴリズムにより実現できる。前方向の投影画像については解析対象としないが、自動設定される眼軸が適切か否かを読影医が判断するための参考にする。

続いて、後述する図７Ａ（ｂ）のように、４つの各投影画像ごとに設定された眼軸上に後極エッジ点Pbより所定の距離の位置に計測用重心Poを定義する。この計測用重心Poは、正常眼の場合は、眼軸Pf-Pbの中点あるいは眼球重心Pgに近くなるが、強度近視眼では眼球の後部側に偏る。これら設定された眼軸および計測用重心を基に、図４で示される４つのパラメータを以下のように解析する。

次いで、ステップＳ３１−１において、上下または左右方向の４枚の投影画像に対し、眼軸を設定し、後極の形状特徴を画像計測する（４投影画像分、繰り返すが、少なくとも１枚の投影画像があれば処理は可能である）。次いで、ステップＳ３１−２において、４枚の投影画像の計測結果を基に、左右方向の対称性、上下方向の対称性、後極尖鋭度を算出する。ここで、左右方向の対称性は上下方向の投影画像の計測結果の平均値を算出し、上下方向の対称性は左右方向の投影画像の計測結果の平均値を算出し、後極尖鋭度は４枚の投影画像の計測結果の平均値を算出する。次いで、ステップＳ５１において、後方向の投影画像に対し、突出数を計測する。
上記の処理により、図４に示したパラメータを実際の眼球に関して求めることができる。

図７Ａ、Ｂは、強度近視の形態解析パラメータの具体例を示す図である。図７Ａ・Ｂをもとに、各パラメータの解析方法の概略について説明する。

（１）眼軸長： 眼軸長は図７Ａ（ａ）をもとに説明する。上方向・下方向・左方向・右方向の４方向の各投影画像ごとに設定された眼軸Ｐ_ｆ−Ｐ_ｂの長さを個別に算出する。続いて、左方向または右方向の投影画像の眼軸Ｐ_ｆ−Ｐ_ｂとＸ軸（水平線）とのなす角αを算出するとともに、上方向または下方向の投影画像の眼軸Ｐ_ｆ−Ｐ_ｂとＹ軸（垂直線）とのなす角βを算出する。左方向または右方向の投影画像から算出される眼軸長に対しては、1/cosβを乗算して上下方向の傾き補正を行い、上方向または下方向の投影画像から算出される眼軸長に対しては、1/cosαを乗算して左右方向の傾き補正を行う。ただし、これらの補正係数による補正量はかなり微小であるため、実用上は省略しても影響ない。最後に、これら４方向の投影画像を基に算出された眼軸長Ｌ_１、Ｌ_２、Ｌ_３、Ｌ_４の平均値（（Ｌ_１＋Ｌ_２＋Ｌ_３＋Ｌ_４）／４）が、所定のしきい値（例．175画素）を超えるか否かで０（正常）または１（異常）と判定する。

（２）後極尖鋭度： 後極尖鋭度は図７Ａ（ｂ）をもとに説明する。上方向・下方向・左方向・右方向の４方向の各投影画像ごとに後極エッジ点Ｐ_ｂにおける曲率を算出する。ただし、曲率半径では固体差が生じにくいので、計測用重心Ｐ_０から後極エッジ点Ｐ_ｂを中心に左右または上下に所定の角度（例えば、±θ_ｎ度）方向に延長した直線の後極輪郭線上の交点をＰ_＋ｎとＰ_−ｎとし、２つの交点とPbとの３頂点のなす角Ｐ_＋ｎ−Ｐ_ｂ−Ｐ_−ｎを算出し、この角度（360度単位）を後極尖鋭度と定義する。具体的には、

により算出できる。そして、４方向の画像を基に算出された後極尖鋭度θ_１、θ_２、θ_３、θ_４の平均値（（θ_１＋θ_２＋θ_３＋θ_４）／４）が、２つの下限しきい値（例．140度）と上限しきい値（例．150度）を設定し、両しきい値範囲であれば０（中間型）、下限しきい値未満は１（紡錘型）、上限しきい値を超えたら２（樽型）と判定する。

（３）左右対称性： 左右対称性は、図７Ｂ（ｃ）を用いて説明する。上方向または下方向の投影画像の各計測用重心Ｐ_０から後極エッジ点Ｐ_ｂを中心に左右に大小２種類の角度（例えば、±θ_ｗ度および±θ_ｎ度）方向に延長した直線の後極輪郭線上の交点をＰ_＋ｗとＰ_−ｗおよびＰ_＋ｎとＰ_−ｎとする。そして、Ｐ_０−Ｐ_ｂ−Ｐ_＋ｗで囲まれた円弧の面積をＳ_＋ｗ、Ｐ_０−Ｐ_ｂ−Ｐ_−ｗで囲まれた円弧の面積をＳ_−ｗ、Ｐ_０−Ｐ_ｂ−Ｐ_＋ｎで囲まれた円弧の面積をＳ_＋ｎ、Ｐ_０−Ｐ_ｂ−Ｐ_−ｎで囲まれた円弧の面積をS_-nとして、４領域の円弧面積を算出する。

続いて、右眼の場合はＳ_＋ｗ×100／Ｓ_−ｗを広領域Ｅ／Ｎ面積比率およびＳ_＋ｎ×100／Ｓ_−ｎを狭領域Ｅ／Ｎ面積比率とし、左眼の場合はＳ_−ｗ×100／Ｓ_＋ｗを広領域E/N面積比率およびＳ_−ｎ×100／Ｓ_＋ｎを狭領域E/N面積比率として算出する。そして、２種類の広領域E/N面積比率と狭領域E/N面積比率より100との差が大きい方のE/N面積比率を左右対称性の判定に採用し、上方向と下方向の投影画像を基に算出された２投影画像のE/N面積比率（耳側／鼻側比率）［％］の平均値を左右対称性とする。最後に、２つの下限しきい値（例．９０％）と上限しきい値（例．１１０％）を設定し、両しきい値範囲であれば０（対称）、下限しきい値未満は１（鼻側非対称）、上限しきい値を超えたら２（耳側非対称）と判定する。

尚、大小２種類の角度範囲でE/N面積比率を算出する理由は、後極部の拡張容積が比較的小さい場合には、後極エッジ点Pbがわずかに左右（通常は耳側）に偏移し、後極部の拡張容積が比較的大きい場合には、後極部が大きく左右（通常は鼻側）に偏移するという強度近視の病態進行上の特性によるものである。前者の場合は狭領域E/N面積比率で非対称性と判定されるが、広領域E/N面積比率では対称と判定されやすい。一方、後者の場合は狭領域E/N面積比率で非対称性（通常は耳側）と判定されるとともに、広領域E/N面積比率でも反対方向（通常は鼻側）に非対称と判定され、後者の非対称性の方が顕著であることを、見出した。

（４）上下対称性：上下対称性は、図７Ｂ（ｄ）を用いて説明する。左方向または右方向の投影画像の各計測用重心Poから、上記左右対称性の算出手法と同様に、後極エッジ点Pbを中心に上下に大小２種類の角度（例えば、±θ_ｗ度および±θ_ｎ度）方向に延長した直線の後極輪郭線上の交点をＰ_＋ｗとＰ_−ｗおよびＰ_＋ｎとＰ_−ｎとする。そして、Ｐ_０−Ｐ_ｂ−Ｐ_＋ｗで囲まれた円弧の面積をＳ_＋ｗ、Ｐ_０−Ｐ_ｂ−Ｐ_−ｗで囲まれた円弧の面積をＳ_−ｗ、Ｐ_０−Ｐ_ｂ−Ｐ_＋ｎで囲まれた円弧の面積をＳ_＋ｎ、Ｐ_０−Ｐ_ｂ−Ｐ_−ｎで囲まれた円弧の面積をＳ_−ｎとして、４領域の円弧面積を算出する。続いて、Ｓ_＋ｗ×100／Ｓ_−ｗを広領域T/B面積比率およびＳ_＋ｎ×100／Ｓ_−ｎを狭領域T/B面積比率として算出する。そして、２種類の広領域T/B面積比率と狭領域T/B面積比率より１００との差が大きい方のT/B面積比率を上下対称性の判定に採用し、左方向と右方向の投影画像を基に算出された２つの投影画像のT/B面積比率（上側／下側比率）［％］の平均値を上下対称性とする。最後に、２つの下限しきい値（例．９０％）と上限しきい値（例．１１０％）を設定し、両しきい値範囲であれば０（対称）、下限しきい値未満は１（下側非対称）、上限しきい値を超えたら２（上側非対称）と判定する。ただし、強度近視では上側非対称と判定される例は稀である。

以上の解析を行うと、眼軸長は０（正常）または１（異常）の２通り、後極尖鋭度は０（中間型）、１（紡錘型）、２（樽型）の３通り、左右対称性は０（対称）、１（鼻側非対称）または２（耳側非対称）の３通り、上下対称性は０（対称）、１（下側非対称）、２（上側非対称）の３通りに分類され、強度近視の病態が理論上は正常眼を含めて全５４通りにタイプ分類されることになる。

以下に、前述した形態解析パラメータの算出を行うための、投影画像に対する具体的な前処理の内容を説明する。

図８（ａ）、図９（ａ）、図１０（ａ）は、ステップＳ３１−１の、上下方向又は左右方向の投影画像の解析処理を示す図である。図８（ｂ）、図９（ｃ）、図１０（ｄ）は、それぞれの処理に対応する画像の概要を示す図である。

まず、ステップＳ３２において、左右方向、上下方向のいずれかの投影画像に対し、投影画像の二値化処理と輪郭抽出する（図８（ｂ）参照）。次いで、ステップＳ３３において、眼軸Ｃの自動設定処理を行う（ステップＳ３３）。ステップＳ３４において、設定された眼軸は適正か否かを判定し、Ｙｅｓの場合には、ステップＳ３６に、Ｎｏの場合にはステップＳ３５に進む。ステップＳ３５では、眼軸の修正処理を行う（図８（ｄ）参照）。ここでは、眼軸Ｃと眼球の輪郭との交点を上下させて眼軸をＣ’に修正する。
図９（ａ）に示すように、ステップＳ３６で、決定された眼軸Ｃで眼軸長（前極側エッジ点から後極側エッジ点までの距離Ｄａ）を計測する（図９（ｂ）： Ｄａ＝｛（Ｘ２−Ｘ１）^２＋（Ｙ２−Ｙ１）^２｝^１／２）。

ステップＳ３７で、眼軸Ｃ上に後極側エッジ点より所定の距離位置（正常眼軸長／２）に計測用の重心点Ｐｏ（Ｘ０，Ｙ０）を追加する。ここで、Ｘ０＝（Ｘ１−Ｘ２）Ｄｓ／Ｄａ＋Ｘ２、Ｙ０＝（Ｙ１−Ｙ２）Ｄｓ／Ｄａ＋Ｙ２である（図９（ｃ）参照）。実施例では、Ｄｓ＝８７ｐｉｘｅｌ（１３．４ｍｍ）である。

ステップＳ３８において、後極側輪郭に、後極側エッジ点を中心に計測用重心点Ｐｏから正負θ_ｗ・θ_ｎ度の方向にエッジ点Ｐ_＋ｗ、Ｐ_＋ｎ、Ｐ_−ｗ、Ｐ_−ｎを追加する（図９（ｄ）参照）。

次いで、図１０に示すように、ステップＳ３９において、後極部−θ_ｎ（例えば２２．５度を設定）〜０度および０度〜＋θ_ｎ度範囲で円弧面積比の算出（狭領域対称性）する（図１０（ｂ）参照）。ステップＳ４０で、後極部−θ_ｗ（例えば４５度を設定）〜０度および０度〜＋θ_ｗ範囲で円弧面積比の算出（広領域対称性）する（図１０（ｃ）参照）。ステップＳ４１で、狭領域面積比と広領域面積比のいずれか比率が大きい方を対称性パラメータとする。次いで、ステップＳ４２において、後極部輪郭の−θ_ｎ度、０度、＋θ_ｎ度の３つのエッジ点を結ぶなす角を算出する（図１０（ｄ）参照： 尖鋭度θ）。

図１１（ａ）は、図５のステップＳ５１からの、後方向投影画像処理の流れを示すフローチャート図である。まず、ステップＳ５２において、後方向（左右各１枚）の投影画像に対し、平滑化フィルタ処理（ノイズ除去）を行う（図１１（ｂ）参照）。次いで、ステップＳ５３において、更新された投影画像に対し、濃淡ピーク点Ｐｎ（ｎ＝１、２、…）をまず１つ抽出する。図１１（ｃ）はピークＰ１を抽出した様子を示す図である。図１１（ｄ）は、ピークＰ２を抽出した様子を示す図である。尚、ピークがない場合には、図１２（ａ）のステップＳ５５に進む。ピークがある場合には、ステップＳ５４に進み、抽出したピーク点Ｐｎの周辺画素値を０にクリアする（図１１（ｅ）、（ｆ）参照）。

図１２（ａ）は、図１１（ａ）に続く処理を示す図である。ステップＳ５４からステップＳ５５に進み、原画像に対し、重心点Ｐｇ（Ｘｇ、Ｙｇ）を算出する（図１２（ｂ）、図１２（ｃ）より、図１２（ｄ）を求める）。

各ピーク点（Ｐｘ１、Ｐｙ１），（Ｐｘ２，Ｐｙ２）の重心Ｐｇ（Ｘｇ，Ｙｇ）からのオフセットと角度（γ１、γ２）を算出する。

次に、レンダリング画像の前処理について図１３、図１４を参照しながら説明する。
図１３（ｂ）は左方向レンダリング画像の例で、画像サイズはW300×H300画素程度（症例によりフレームからはみ出すこともあるため、画像サイズは若干増減する）で、階調はモノクロ256階調である。これに対して、固定しきい値（例えば、0〜255値に対して25を設定）で二値化を行った結果が図１３（ｃ）であり、眼球内部は１、背景は０で与えられる。この時、原画像の撮像状態により階調が均一でないと、右下部のように、二値化によりゴミ領域が残留する場合がある。数画素のゴミについては、上記特許文献８に記載のように、二値化処理の前に３×３の平均化マトリックス等を用いて、事前に平滑化を加えることで回避できる。比較的面積の大きいゴミ領域に対しては、二値画像から８近傍画素の連結性に基づいて閉領域を抽出する処理を行い、面積が所定のしきい値以下（例えば、200画素に設定）の閉領域を削除する処理を行う。その結果、眼球内部は２、背景は０で与えられる。最後に、眼球内部と背景との境界部分に輪郭画素を設定する。具体的には眼球内部に設定されている画素の８近傍画素のうち、４画素以上が背景画素になっている場合、輪郭画素として画素の値を１に変更する。その結果が図１３（ｄ）であり、眼球内部は２、輪郭線上は１、背景は０で与えられる。

より詳細に説明すると、図１３は、投影画像の二値化処理と輪郭抽出処理（Ｓ３２）の詳細について説明するための図である。まず、ステップS３２−１において、左右方向、上下方向のいずれかの投影濃淡画像に対し、平滑化フィルタ処理を行う（図１３（ｂ）参照）。
ここで、ソース画像データ２５６階調であり、Image（x，y）＝０〜２５５ (x=0,…,Xs-1; y=0,…,Ys-1)である。

ステップS３２−２において、投影画像の二値化処理を行う（図１３（ｂ）参照）。
ここで、二値画像データは、
ImageM(x,y)={0：背景，1：内部}である。

次いで、ステップＳ３２−３において、近傍画素の連結処理、領域抽出処理を行う。
ここで、多値画像データは、
ImageM(x,y)={0：背景，2以上：内部} (x=0,…,Xs-1; y=0,…,Ys-1) である。

次いで、図１４（ａ）のステップＳ３２−４において、近傍領域の統合処理を行う（図１４（ｂ）参照）。
ここで、多値画像データは、
ImageM(x,y)={0：背景，2以上：内部} (x=0,…,Xs-1; y=0,…,Ys-1)である。

次いで、ステップＳ３２−５において、微小面積領域の削除処理を行う（図１４（ｃ）参照、３の領域が削除される）。
ここで、多値画像データは、
ImageM(x,y)={0：背景，2：内部} (x=0,…,Xs-1; y=0,…,Ys-1) である。

次いで、ステップＳ３２−６において、輪郭抽出処理を行う（図１４（ｄ）参照）。
ここで、三値画像データは、
Image(x,y)={0：背景，1:輪郭，2：内部} (x=0,…,Xs-1; y=0,…,Ys-1) である。このようにして、眼球の投影画像を、背景、輪郭、内部の３値に区別できる。

図１５Ａは、眼軸の自動設定アルゴリズムの概要を示す図である。
前述の通り、眼軸は後極の左右対称性や上下対称性などのパラメータを算出する際の基準になるもので、高精度に設定されることが要求される。しかし、病態により眼球の形状自体が大幅に変化する以上、画一的な眼球形状特徴を手がかりにして自動的に眼軸を設定することは困難である。強度近視では一般に後極側が非対称になっても、前極側は対称性が維持されるということを仮定し、前極の対称性を基に眼軸を自動設定する手法を提案する。そして、前極も非対称である場合など、この仮定が成立しない場合には、手動による対話型修正を加えて個別に対応するものとする。

図１４（ｄ）の画像に対して、値が１となる画素のＸ方向の最小値Ｘ_１と最大値Ｘ_２を求め、各々最小最大をとるときのＹ方向の画素位置をＹ_１およびＹ_２とする。また、値が１以上となる画素のＸ座標の平均値をＸｇ、Ｙ座標の平均値をＹｇとする。これらの座標値を基に、図１５Ａ（１）に示されるように、(Ｘ_１,Ｙ_１)を前極エッジ点Ｐｆ、(Ｘ_２,Ｙ_２)を後極エッジ点Ｐｂ、(Ｘｇ,Ｙｇ)を重心点Ｐｇと定義する。この時、Ｐｆ−Ｐｇ−Ｐｂが眼軸に対応するが、正常眼のように前極も後極も対称性をもつ眼球でない限り、この段階で仮設定される眼軸は適切な方向にならない。

Ｐｆ, Ｐｇ, Ｐｂの３点で最も信頼できる点はＰｇである。そこで、Ｐｇを基準にしてＰｆを補正する。具体的には、図１５Ａ（２）のように、前極エッジ点ＰｆをＰｆのＸ座標がＰｇに近づかない範囲で、前極の輪郭線上に上下に微少量移動させる。移動させた各位置で、Ｐｇ−Ｐｆの直線上にＰｃなる中間計測点を設定し、Ｐｃより直線Ｐｇ−Ｐｆと直交する上下方向に直線を延ばし、輪郭線との交点をＰ_＋ｖおよびＰ_−ｖとし、距離Ｐｃ−Ｐ_＋ｖとＰｃ−Ｐ_−ｖとの差を計算する。この差が最小値をとるときのＰｆの位置を前極エッジ点に決定する。このとき、Ｐｃは一箇所でなく、直線Ｐｇ−Ｐｆ上に所定の範囲で動かし（例えば、ＰｃをＰｇからＰｆ方向に２〜１２画素の距離だけ動かす）、各々算出される距離Ｐｃ−Ｐ_＋ｖとＰｃ−Ｐ_−ｖとの差の総和が最小値をとるＰｆの位置を探索する。続いて、逆に決定されたＰｆを基準にしてＰｇを補正する。具体的には、図１５Ａ（３）のように、重心点ＰｇをＰｇのＹ座標が上下の輪郭線に近づかない範囲で、Ｘ座標Ｘｇを固定にして上下に微少量移動させる。上記と同様に、移動させた各位置で、Ｐｆ−Ｐｇの直線上にＰｃなる中間計測点を設定し、Ｐｃより直線Ｐｆ−Ｐｇと直交する上下方向に直線を延ばし、輪郭線との交点をＰ_＋ｖおよびＰ_−ｖとし、距離Ｐｃ−Ｐ_＋ｖとＰｃ−Ｐ_−ｖとの差を計算する。この差が最小値をとるときのＰｇの位置を重心点に決定する。同様に、Ｐｃは一箇所でなく、直線Ｐｆ−Ｐｇ上に所定の範囲で動かし（例えば、ＰｃをＰｆからＰｇ方向に２〜１２画素の距離だけ動かす）、各々算出される距離Ｐｃ−Ｐ_＋ｖとＰｃ−Ｐ_−ｖとの差の総和が最小値をとるＰｇの位置を探索する。

最後に、決定されたＰｆとＰｇを基準にして、Ｐｆ−Ｐｇ−Ｐｂが一直線上に近づくようにＰｂを補正する。具体的には、図１５Ａ（４）のように、後極側エッジ点ＰｂをＰｂのＸ座標がＰｇに近づかない範囲で、後極の輪郭線上に上下に微少量移動させる。移動させた各位置で、重心点Ｐｇから直線Ｐｆ−Ｐｂに下ろした垂線と直線Ｐｆ−Ｐｂとの交点と重心点Ｐｇとの距離を算出する。この距離が最小値をとるときのＰｂの位置を後極エッジ点に決定する。

以上の処理により、Ｐｆ−Ｐｂを眼軸と決定できるが、眼軸位置を画面に表示して読影医による判断を仰ぎ、不適切な場合、２点PfとPbを画面上で対話形式に修正し、修正されたＰｆとＰｂの中点をＰｇに再設定する。また、パラメータ解析を行うため、図７Ａ（ｂ）に示されるように、決定された眼軸上でPbより所定の距離（例えば、87画素）位置に計測原点Poを設定し、 Poから後極エッジ点Pbを中心に上下に大小２種類の角度（例えば、±θ_ｗおよび±θ_ｎ）方向に延長した直線の後極輪郭線上の交点Ｐ_＋ｗとＰ_−ｗおよびＰ_＋ｎとＰ_−ｎを設定する。

より詳細な処理内容を以下に説明する。
図１５Ｂは、眼軸の自動設定処理（ステップＳ３３）の詳細を示す図である。
図１５Ｂ（ａ１）に示すように、ステップＳ３３−１において、輪郭抽出後の投影画像に対し、前極側および後極側のエッジ点（Ｘ１，Ｙ１）、（Ｘ２，Ｙ２）と重心点（Ｘｇ，Ｙｇ）を算出する。左右方向の投影画像の場合（図１５Ｂ（ａ２））、Ｘ軸方向の最大最小エッジ点を前後極側のエッジ点とし、上下方向の投影画像の場合、Ｙ軸方向の最大最小エッジ点を前後極側のエッジ点とする（ステップＳ３３−２）。図１５Ｂ（ｂ）から図１５Ｂ（ｅ）までは、左右、上下方向の投影画像に関する処理を行った様子を示す図である。

図１６（ａ）は、図１５に続く処理の流れを示す図である。ここで、前極側は対称性が保持されていることを前提とする。ステップＳ３３−３において、前極側エッジ点の補正と決定処理を行う（図１６（ｂ）参照）。ここで、前極側エッジ点を輪郭線上で微少量動かし、重心点Ｐｇから前極側エッジ点方向のベクトルから垂直方向に上下に輪郭線まで延ばした２つの距離（両矢印）が最も近くなるように、前極側エッジ点Ｐｆを決定する。

次いで、ステップＳ３３−４において、重心点の補正と決定処理を行う（図１６（ｃ）参照）。
重心点Ｐｇを輪郭線内部で垂直方向に上下に動かし、既に決定されている前極側エッジ点Ｐｆから重心点方向のベクトルから垂直方向に上下に輪郭線まで延ばした２つの距離（両矢印）が最も近くなるように、重心点Ｐｇ’を決定する。

次に、図１７（ａ）に示すように、ステップＳ３３−５において、後極側エッジ点の補正と眼軸の決定処理を行う（図１７（ｂ）参照）。ここで、後極側エッジ点を輪郭線上で微少量動かし、既決定の前極側エッジ点Ｐｆから重心点Ｐｇ’のベクトルの延長線上に最も近くなるように、後極側エッジ点Ｐｂを決定する。

図１８（ａ）は、眼軸の修正処理の流れを示す図である。まず、ステップＳ３５−１において、画面上で前極側エッジ点Ｐｆまたは後極側エッジ点Ｐｂのいずれかを選択する（図１８（ｂ）参照）。

次いで、ステップＳ３５−２において、選択した前極側エッジ点Ｐｆまたは後極側エッジ点Ｐｂを画面上で所望の位置（Ｐｘ）へ移動させる（図１８（ｃ）参照）。

次いで、移動させて設定したエッジ点Ｐｘの輪郭線ＬＲ上への自動補正と重心点Ｐｇ’の自動補正を行う（図１８（ｄ）参照）。エッジ点を最も近傍の輪郭線上の位置へ移動させ、更新したエッジ点ともう一方の更新しないエッジ点を結ぶ直線上に重心点を移動させる）。

ステップＳ３５−４において、適切な方向か否かを判定し、ＹＥＳであれば、処理を終了し、Ｎｏであれば、ステップＳ３５−１に戻る。

図１９（ａ）は、後極輪郭部における計測用エッジ点の作成処理（Ｓ３８）の詳細を示す図である。
ステップＳ３８−１において、計測用重心点（Ｘ０，Ｙ０）より眼軸Ｃと直交する方向に上下の輪郭線ＬＲ上に９０度および−９０度計測用エッジ点（Ｘｐ９，Ｙｐ９）、（Ｘｍ９，Ｙｍ９）を追加する（図１９（ｂ）参照）。

ステップＳ３８−２において、９０度、０度、−９０度計測用エッジ点をもとに、θ_ｗおよび−θ_ｗ度計測用エッジ点（Ｘｐ４，Ｙｐ４）、（Ｘｍ４，Ｙｍ４）を追加する（図１９（ｃ）参照）。

ステップＳ３８−３において、θ_ｗ、０度、−θ_ｗ計測用エッジ点をもとに、θ_ｎおよび−θ_ｎ度計測用エッジ点（Ｘｐ２，Ｙｐ２）、（Ｘｍ２，Ｙｍ２）を追加（図１９（ｄ）参照）。
このようにして後極部に計測用エッジ点を作成する。

図２０（ａ）は、平滑化フィルタ処理（図１１（ａ）：Ｓ５２）の流れを示すフローチャート図である。図２０（ｂ）は、平滑化フィルタの例を示す図である。

フィルタ演算は、マトリックスを画像データに対して掛け合わせて処理するものである。処理する注目する画素の近傍を領域を複数個選んで参照する（例えば３×３、上下左右）。この値（１／９）を画素の値に乗算して、加算する処理により平均値を取る（３×３の９画素の総和の平均値をとる）。

図２０（ｂ）に示すように、３×３の平滑化フィルタ・マトリックスを、M（ｉ，ｊ）の例として示す図である。ステップＳ５２−１において、図２０（ｂ）に示す３×３の平滑化マトリックスを定義する。次いで、ステップＳ５２−２において、ｘ＝ｙ＝１，ｉ＝ｊ＝−１、一時的な変数（画素値の一時的な変数）ｖ＝０として初期化の定義とを行う。ｖはｘ，ｙは、画素をスキャンするｘ方向のアドレスとｙ方向のアドレス（トータルＸｓ、Ｙｓ個の画素）の先頭を決める。ｉ、ｊは、マトリックスを演算するために加える変位であり、入力の投影画像は、Ｉｍａｇｅ（ｘ、ｙ）であり、ｘ＝（０，１，…，Ｘｓ−１）、ｙ＝（０，１，…，Ｙｓ−１）で定義できる。ｉとｊとは、−１から＋１の範囲の整数値を取る。

ステップＳ５２−３で、v←v+Image(x+i,y+j)・M(i,j)とする。ステップＳ５２−４で、ｉを１ずつインクリメントし、ステップＳ５２−５で、ｉと＋１とを比較し、ｉ＞１になれば、１つの３×３のフィルタによる処理を終えて、ステップＳ５２−６に進み、ｊに＋１を行い、同様の処理を行うことで、３×３のマトリックスを乗算する処理を終える。そして、ステップＳ５２−９、ステップＳ５２−１１で、全ての画素に関する演算が行われたか否かを判定し、処理を終了する（Ｅｎｄ）。このようにして、３×３の平滑化マトリックスを入力投影画像上で動かすことにより、３×３の平滑化フィルタにより、投影画像の画素値の平滑化を行うことができる。

図２１は、投影画像に対する濃淡ピーク点の抽出処理（図１１（ａ）：Ｓ５３参照）の流れを示すフローチャート図である。ｘとｙとでループ処理を行う点は、図２０と同様である。但し、Ｘｓ、Ｙｓまでの処理を行う。画素Ｉｍａｇｅ（ｘ，ｙ）は、０〜２５５までの値を持っている。投影画像の中で局所的にピークをもつローカルピークを見い出す処理である。画素が０から２５５までの値をもっており、ステップＳ５３−１において、初期値の例として１２８より小さくない値としてＶｍ＝１２８に設定し、ピーク座標（Ｐｘ，Ｐｙ）として、初期値（−１，−１）を、設定する。また、濃淡画像Ｉｍａｇｅ（ｘ，ｙ）をコピーした一時的なコピーＩｍａｇｅ２（ｘ，ｙ）をバックアップとして持つとともに、初期値としてｘ＝ｙ＝０を設定する。ステップＳ５３−２において、Ｉｍａｇｅ２（ｘ，ｙ）と濃淡のピークＶｍとを比較する。Ｉｍａｇｅ２（ｘ，ｙ）＞Ｖｍであれば、ステップＳ５３−３に進み、Ｖｍ＝Ｉｍａｇｅ２（ｘ，ｙ）とし、Ｐｘ＝ｘ，Ｐｙ＝ｙを保存する。以下、全ての画素において処理を行うと、ピーク点が求まる。すなわち、Ｉｍａｇｅ２（ｘ，ｙ）＜＝Ｖｍであれば、ステップＳ５３−３をスルーし、ステップＳ５３−４で、ｘを１だけインクリメントする。画像の水平方向の末端にあたるまで（ｘ＝Ｘｓ）処理を繰り返し（ステップＳ５３−５）、次いで、ｙについても処理を行う（ステップＳ５３−６、７）。ステップＳ５３−８で、Ｘｐ＜０かつＹｐ＜０であるか否かを判定し、Ｙｅｓであれば、濃淡のピークが存在しないとして（ピーク点無し）、例えば初期値を再設定する。Ｎｏであれば、ピークあり（Ｘｐ，Ｙｐ）、このピーク（Ｘｐ，Ｙｐ）を出力する。

次いで、図２２に進む。図２２は、図２１の処理で抽出できたピーク点の周辺画素のクリア処理の流れを示すフローチャート図である（図１１のステップＳ５４）。図２１のピーク出力から、ステップＳ５４−１において、ピーク座標（Ｘｐ，Ｙｐ）で、ｘ＝ｙ＝０としてリセットする。ステップＳ５４−２において、ピーク点（Ｘｐ，Ｙｐ）までの距離ｒを
ｒ＝{(x-Px)²+(y-Py)²}^1/2
により算出する。ステップＳ５４−３において、距離ｒと、そのしきい値（例えば５０画素）とを比較する。ｒがＳｒ以下であれば、ステップＳ５４−４に進み、Image2(x,y)=0として、ｒがＳｒより大きい場合とともに、ステップＳ５４−５でｘを１だけインクリメントし、ステップＳ５４−６において、ｘとＸｓとを比較する。ｘ＜Ｘｓであれば、ステップＳ５４−２に戻る。ｘがＸｓ以上であれば、ステップＳ５４−７に進み、ｙについても同様の処理を行い、クリア処理を終了する。これにより、ピーク点（Ｘｐ，Ｙｐ）から半径Ｓｒの範囲の領域における画素値をクリア処理することができる。

図２３は、原画像に対する重心点の算出処理を示すフローチャート図である（図１２（ａ）：ステップＳ５５）。ここでは、画素の平均値を取ることで、重心点を求めている。

まず、ステップＳ５５−１で、重心座標(Xg,Yg)において, Xg=Yg=0, 重み総和V=0, x=y=0と初期化する。次いで、ステップＳ５５−２において、
Xg←Xg+Image(x,y)・x, Yg←Yg+Image(x,y) ・y, V←V+ Image(x,y)
として、ステップＳ５５−３からステップＳ５５−６までの処理により、画素の平均値を計算し、ステップＳ５５−７で、
Xg←Xg/V, Yg←Yg/V
と、原画像に対する重心点を求めることができる（図１２（ｄ）参照）。

図２４は、画素の二値化処理（ステップＳ３２）の詳細を示すフローチャート図である。画素値０〜２５５を０か１かのいずれかに変換する。すなわち、Sl（エル）よりも大きければ“１”にして、そうでなければ“０”にする。

まず、ステップＳ３２−２−１において、二値画像データImageM(x,y)=0（初期値は全画素０）,ｘ＝ｙ＝１を設定する。ステップＳ３２−２−２において、Image(x,y)をしきい値Sl（エル）と比較する。しきい値Slは、例えば２５である。Image(x,y)がしきい値Slより大きい場合には、ステップＳ３２−２−３に進み、二値画像データImageM(x,y)=1とし、Image(x,y)がしきい値Sl（エル）以下である場合とともに、ステップＳ３２−２−４に進み、ｘを１だけインクリメントし、ステップＳ３２−２−５に示すように、ｘがＸｓ−１になるまで上記の処理を繰り返す。ｘがXs−１より大きくなると、ステップＳ３２−２−６に進み、ｙについて１だけインクリメントし、ｙがＹｓ−１より大きくなるまで処理を繰り返す。ｙがYs−１以上になると、処理が終了する。これにより、Sl（エル）よりも大きければ“１”にして、そうでなければ“０”にする画像の二値化処理をすることができる（図１３（ｃ）参照）。

図２５（ａ）は、近傍画素の連結処理、領域抽出処理の流れを示すフローチャート図である（ステップＳ３２−３、図１３（ｄ）参照）。図２５（ｂ）は、８近傍画素のうちの４つの近傍画素を示す図である。０、１で２値化した場合に、例えば、対象画素の値が１で、２以上の値（領域Ｉｄ）を持ったものが対象画素の４つの近傍画素のいずれかに存在するかをみる。２以上の値をもつ近傍画素が存在する場合、対象画素の値を当該近傍画素のＩｄに設定するとともに、当該近傍画素のＩｄに対応する領域テーブル内の面積値（領域内の画素の個数）を１だけインクリメントする。２以上の値をもつ近傍画素が存在しない場合、対象画素に新規な領域Ｉｄ（初期値は２）を割り当てるとともに、当該Ｉｄに対応する領域テーブルに新規なレコードを登録し、その初期面積値を１に設定するする。このような処理を繰り返すことにより、０、１で２値化された画像の中で１をもつ全画素に２以上の領域Ｉｄが割り当てられ、Ｎ個の領域に分割されるとともに、領域テーブルに各領域の面積算出結果が格納される。

まず、ステップＳ３２−３−１において、Id: 領域Id (初期値:1), 領域テーブルTb(id)：初期値0, n:領域数(初期値ｎ＝０)、ｘ＝ｙ＝１を設定する。ステップＳ３２−３−２において、対象画素（ｘ、ｙ）が眼球内部にあるか否か、即ち、ImageM(x,y)と０とを比較する。眼球内部にある場合、即ち、０より大きい場合には、ステップＳ３２−３−３において、４つの８近傍画素(x-1,y-1), (x,y-1), (x+1,y-1), (x-1,y)を確認する。すなわち、自分のアドレスよりもアドレスが小さい４つの画素について、背景（０）と眼球内部（２以上、アドレスが小さい既処理画素では値１は存在しない）のうち、内部と指定されているものがあるか否かを確認する。即ち、ステップＳ３２−３−４において、２以上の値Idpをもつ画素があるか否かを判定する。

Ｙｅｓの場合は近傍画素のIdpを自己の画素の領域Ｉｄとして流用し、Ｎｏの場合には、ステップＳ３２−３−６に進み、新規領域Ｉｄの割り当てを行い、Id←Id+1, Idp=Id, n←n+1とし、Ｙｅｓの場合とともに、ステップＳ３２−３−５に進み、テーブルの値（２番）に画素のカウント値（初期値０）に１を加算して１とし、同時に領域番号Ｉｄｐ（＝２）として画像データを更新する。
Tb(Idp)←Tb(Idp)+1, ImageM(x,y)=Idp

次いで、ステップＳ３２−３−７に進み、ｘを１だけインクリメントし、ステップＳ３２−３−８において、ｘとＸｓ−１とを比較し、輪郭まで処理を繰り返し、次いで、ｙについても同様の処理を行う。これにより、近傍画素の連結処理、領域抽出処理を終了する（Ｅｎｄ：図１３（ｄ）参照）。これにより、領域分割と各領域の面積計算が完成する。しかし、領域Ｉｄ割り当て時に、４近傍画素しか参照していないため、本来同一Ｉｄが割り当てられるべき近傍の領域に異なる領域Ｉｄが割り当てられ、過剰に領域分割されがちになる。そこで、近傍領域を統合する補正処理を行う。

図２６（ａ）は、近傍領域の統合処理の流れを示すフローチャート図である（図１４、ステップＳ３２−４、図１４（ｂ）参照）。図２６（ｂ）は、８近傍画素の例を示す図である。領域テーブルが完成した状態で、８近傍画素全部を見て抽出した領域を整理する処理であり、全画素の領域番号が再割り当てされて整理される。

上記処理を行った結果、多値画像データは、
ImageM(x,y)={０：背景，２以上：内部} (x=0,…,Xs; y=0,…,Ys)
になっており、２以上の値をもつ画素をできるだけ小さい値（２以上の範囲で）になるように統合する。

まず、ステップＳ３２−４−１で、ｘ＝ｙ＝１とし、ステップＳ３２−４−２において、対象画素（ｘ、ｙ）が領域内部にあるか否か、即ち、ImageM(x,y)が正か０かを判定する。正の場合には、ステップＳ３２−４−３に進み、８近傍画素(x-1,y-1), (x,y-1), (x+1,y-1), (x-1,y), (x+1,y), (x-1,y+1), (x,y+1), (x+1,y+1)を確認する。すなわち、ステップＳ３２−４−４において、ImageM(x+i,y+j)< ImageM(x,y)となる近傍の異なる領域内画素が存在するか否かを判定する。即ち、２以上の値が入っていて、自分のＩＤ番号よりも若い番号の付けられている近傍の画素を調べる。Ｙｅｓの場合には、ステップＳ３２−４−５に進み、領域Ｉｄ更新処理を行い、Ｎｏの場合、ステップＳ３２−４−２でImageM(x,y)が０の場合と同様に、ステップＳ３２−４−６に進み、ｘを１だけインクリメントし、処理を継続する。ｘがＸｓ−１になると、ｙに関して同様の処理を行い、ｙ＝Ｙｓ−１になると、処理を終了する（Ｅｎｄ）。

次いで、図２７に示すように、自分のＩＤ番号よりも若い番号の付けられている近傍の画素を調べて若い番号が付いている場合には、自分の画素を若い番号の方に入れ直す（若いグループに入れる）。自分の領域ｉｄを周辺画素(x+i,y+j)のＩｄｐに入れ替える場合、画像データ全体に渡って、自分と同一の領域ｉｄをもつ全画素をＩｄｐに変更する。そうすると、当該領域ｉｄの値をもつ画素が無くなるため、当該領域ｉｄの値より大きい値をもつ全ての画素の値を１だけ減少させる。これにより、使用される領域ｉｄの種類Ｎが１だけ減少し、本処理を繰り返すことにより、領域が順次統合されＮが小さくなる。

領域ステップＳ３２−４−１１において、
ｉd= ImageM(x,y), Idp=ImageM(x+i,y+j), u=v=1を設定し、ステップＳ３２−４−１２において、ImageM(u,v)とｉdとを比較する。ImageM(u,v)が０の場合には、ステップＳ３２−４−１５に進む。ImageM(u,v)がｉdの場合には、ステップＳ３２−４−１３に進み、ImageM(u,v)=Idpとし、ステップＳ３２−４−１５に進む。ImageM(u,v)＞ｉdの場合には、ステップＳ３２−４−１４に進み、ImageM(u,v)←ImageM(u,v)-1とし、ステップＳ３２−４−１５に進む。ステップＳ３２−４−１５においては、ｕを１だけインクリメントし、ステップＳ３２−４−１５においてｕとＸｓ−１とを比較し、u≧Xs-1であれば、ステップＳ３２−４−１７に進み、u<Xs-1であればステップＳ３２−４−１２に戻る。ステップＳ３２−４−１７においては、u=1, v←v+1として、ステップＳ３２−４−１８において、ｖとＹｓ−１とを比較する。v≧Ys-1の場合には、ステップＳ３２−４−１９に進み、v≧Ys-1の場合にはステップＳ３２−４−１２に戻る。ステップＳ３２−４−１９においては、領域テーブルTbに関して、領域ｉｄをＩｄｐに統合させ、ｉｄを０クリアする。即ち、
Tb(Idp)←Tb(Idp)+Tb(ｉd), Tb(ｉd)=0, id←ｉd+１
とし、ステップＳ３２−４−２０において、idとn+2とを比較し、id≧n+2であれば、ステップＳ３２−４−２２に進み、id＜n+2であれば、ステップＳ３２−４−２１に進む。ステップＳ３２−４−２１では、全ての領域のｉｄ番号を１だけ減少させる。即ち、
Tb(id-1)←Tb(id), id ← id+1
として、ステップＳ３２−４−２０に戻る。最終的にid≧n+2になると、ステップＳ３２−４−２２において、ｎを１だけデクリメントして処理を終了する。これにより図１４（ｂ）に示すような、近傍領域の統合処理を行うことができる。

図２８は、周辺のごみのような微小面積を削除して、眼軸の決定における間違いを低減するための、微小面積領域の削除処理（ステップＳ３２−５）を示すフローチャート図である（図１４（ｃ）参照）。まず、ステップＳ３２−５−１においてid＝２とし、ステップＳ３２−５−２において、Tb(id)とStとを比較する。ここで、Stは、面積の下限しきい値であり、例えば２００である。Tb(id)＞＝Stであれば、ステップＳ３２−５−１０に進む。Tb(id)＜Stであれば、ステップＳ３２−５−３に進み、画像データ中の全てのidという値をもつ画素を背景画素に設定する。ステップＳ３２−５−３において、変数u、vを１とし、ステップＳ３２−５−４において、ImageM(u,v)と領域番号idとを比較する。ImageM(u,v)＝idであれば、ステップＳ３２−５−５において、ImageM(u,v)=0として領域内の画素を背景領域にしていき、ステップＳ３２−５−６に進む。ImageM(u,v)がidと等しくなければ何もせずに、ステップＳ３２−５−６に進む。

ステップ３２−５−６において、uを１だけインクリメントし、ステップ３２−５−７において、uとXs-1とを比較する。u＜Xs-1であればステップＳ３２−５−４に戻り、u＞＝Xs-1であればステップＳ３２−５−８において、u=1, v←v+1とする。

次いで、ステップＳ３２−５−９において、vとYs-1とを比較し、v≧Ys-1であればステップＳ３２−５−４に戻り、v≧Ys-1であればステップＳ３２−５−１０に進み、Idを１だけインクリメントし、ステップＳ３２−５−１１において、idとn+2とを比較する。id<n+2であれば、ステップＳ３２−５−４に戻り、id≧n+2であれば、ステップＳ３２−５−１２に進む。

このように、前の処理で面積の小さい領域を削除した後に、まだ領域番号３以上の値が残った場合に、無条件に３以上の値を２に変更する（３値化する、図１４（ｃ）参照）。すなわち、ステップＳ３２−５−１２においてu=v=1として、ステップＳ３２−５−１３において、ImageM(u,v)と2との大小を比較し、ImageM(u,v)＞2であれば、ステップＳ３２−５−１４に進み、ImageM(u,v)=2とし、ImageM(u,v)＜＝2の場合とともに、ステップＳ３２−５−１５において、ｕを１だけインクリメントする。これにより、図１４（ｃ）の３以上の領域を２に変更する。ステップＳ３２−５−１６において、uとXs-1とを比較しu<Xs-1であれば、ステップＳ３２−５−１３に戻り、u＞＝Xs-1であれば、ステップＳ３２−５−１７に進み、u=1, v←v+1とする。次いで、ステップＳ３２−５−１８において、vとYs-1とを比較し、v＜Ys-1であれば、ステップＳ３２−５−１３に戻り、v＞＝Ys-1であれば、処理を終了する。これにより、あるしきい値Ｓｔよりも小さい領域を削除する（図１４（ｂ）の３番の領域を背景（０）に置き換える処理を行うこと、すなわち微小面積の領域を削除することができる。

図２９は、輪郭抽出処理の流れを示すフローチャート図である（図１４（ｄ）参照）。自分自身も含めて８近傍の合計９つの画素において、２（内部）になっている画素の数を数え、それが３番から７番の値であれば輪郭として、“１”に変更する処理である。ステップＳ３２−６−１において、ｘ＝ｙ＝１を設定する。次いで、画素カウンタにおいて、ｃ＝０、ｉ＝ｊ＝−１を設定する。ステップＳ３２−６−３において、対象画素（ｘ、ｙ）の８近傍画素(x+i,y+j)が領域内部にあるか否か、即ち、ImageM(x+i,y+j)と0とを比較し、ImageM(x+i,y+j)＞0であれば、ステップＳ３２−６−４において、ｃを１だけインクリメントし、ステップＳ３２−６−３において、＝０の場合とともに、ステップＳ３２−６−５において、ｉを１だけインクリメントする。次いで、ｉと１とを比較し、i≦1であれば、ステップＳ３２−６−３に戻る。ｉと１とを比較し、i＞1であれば、ステップＳ３２−６−６に進み、ｉを１だけインクリメントする。ステップＳ３２−６−６において、ｉと＋１とを比較し、i≦1であればステップＳ３２−６−３に戻り、i>1であれば、ステップＳ３２−６−７において、i=-1, j←j+1とする。

次いで、ステップＳ３２−６−８において、jと+1とを比較し、j≦1であれば、ステップＳ３２−６−３に戻り、j>1であれば、ステップ３２−６−９において、カウントした領域内部にある近傍画素の個数ｃの値を基に対象画素が輪郭になるかを判定する。例えば、3≦c≦7であるか否か（８近傍画素の数が３から７の範囲に入っているかで輪郭と判定する。ｙｅｓの場合には、ステップＳ３２−６−１０において、ImageM(x,y)=1と輪郭設定し、ｎｏである場合とともに、ステップＳ３２−６−１１において、ｘを１だけインクリメントする。次いで、ステップＳ３２−６−１２において、ｘとＸｓ−１を比較し、x<Xs-1であればステップＳ３２−６−２に戻る。x≧Xs-1であれば、ステップＳ３２−６−１３に進み、x=1, y←y+1とし、ステップＳ３２−６−１４において、yとYs-1とを比較し、y≧Ys-1であればステップＳ３２−６−２に戻る。y≧Ys-1になると、処理を終了することで、輪郭を抽出することができる。

図３０は、前極側および後極側エッジ点と重心点の算出処理の流れを示すフローチャート図である（Ｓ３３−１：図１５Ｂ）。前述する重心を求める処理と同様であるが、０から２５５までの階調画像ではなく、０、１、２の３つの値のみをとる。

まず、ステップＳ３３−１−１において、重心座標(Xg,Yg), Xg=Yg=0, 重み総和V=0, x=y=0、前極側エッジ点(X1,Y1), X1=Y1=-1,後極側エッジ点(X2,Y2), X2=Y2=-1を設定する。

次いで、ステップＳ３３−１−２において、対象画素（ｘ、ｙ）が領域内部または輪郭にあるか否か、即ち、ImageM(x,y)と0とを比較する。ImageM(x,y)＝０の場合には、ステップＳ３３−１−１３に進む。ImageM(x,y)＞０の場合には、ステップＳ３３−１−３に進み、Xg←Xg+x, Yg←Yg+y, V←V+1とする。次いで、ステップＳ３３−１−４において、左右方向の投影画像であるか否かを判定する。Ｙｅｓの場合には、ステップＳ３３−１−５に、Ｎｏの場合にはステップＳ３３−１−９に進む。

ステップＳ３３−１−５において、X1<0 or x<X1であるか否かを判定し、Ｙｅｓの場合には、ステップＳ３３−１−６においてX1=x, Y1=y（最小値候補）とし、Ｎｏの場合とともにステップＳ３３−１−７に進み、X>X2であるか否かを判定し、Ｙｅｓの場合にはX2=x, Y2=yとし。Ｎｏの場合と共に、X2=x, Y2=y（最大値候補）とする。このようにして、ｘ方向で最小と最大を見つける処理である。

ステップＳ３３−１−４でＮｏの場合には、上下方向についても同様に、ステップＳ３３−１−９に進み、Y1<0 又は y<Y1であるか否かを判定し、Ｙｅｓの場合にはステップＳ３３−１−１９に進み、X1=x, Y1=yとし、Ｎｏの場合とともに、ステップＳ３３−１−１１において、ｙ＞Ｙ２であるか否かを判定する。Ｙｅｓの場合には、ステップＳ３３−１−１２において、X2=x, Y2=yとし、Ｎｏの場合とともに、ステップＳ３３−１−１３において、ｘに１をインクリメントする。ステップＳ３３−１−１４において、ｘとＸｓとを比較し、x<Xsであれば、ステップＳ３３−１−２に戻り、x≧XsであればステップＳ３３−１−１５においてx=0, y←y+1とし、ステップＳ３３−１−１６において、yとYsとを比較する。y≧Ysであれば、ステップＳ３３−１−２に戻り、y≧Ysであれば、ステップＳ３３−１−１７において、Xg←Xg/V, Yg←Yg/Vとして、ステップＳ３３−１−１８において、右・下方向投影であるか否かを判定し、Ｎｏであれば、そのまま終了し、Ｙｅｓであれば上下左右をひっくり返して、ステップＳ３３−１−１９において、X1≪X2, Y1≪Y2として処理を終了する。これにより、Ｘ１、Ｘ２、Ｙ１、Ｙ２を求めることができる。
仮決定された点に順番に補正する。

図３１は、前極側エッジ点の補正と決定処理の流れを示すフローチャート図である（Ｓ３３−３、図１６参照）。ステップＳ３３−３−１において、Ｘ１、Ｙ１を補正するために、(X1,Y1)と(Xg,Yg)より、後述する方法により距離差dを算出する。まず、ステップＳ３３−３−２において、前極側エッジ点(x1,y1)、初期値x1=X1, y1=Y1、最小距離差dm=d, x=y=0を設定する。

ステップＳ３３−３−３において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、ImageM(x,y)と1とを比較し、異なれば、輪郭上の点でないと判断しステップＳ３３−３−８に進み、同じであれば、輪郭上の点であると判断しステップＳ３３−３−４に進み、
左方向：x<{X1*(Sw-1)+Xg}/Sw
右方向：x>{X1*(Sw-1)+Xg}/Sw
上方向：y<{Y1*(Sw-1)+Yg}/Sw
下方向：y>{Y1*(Sw-1)+Yg}/Sw
を満たすか否かを判定する。ここで、Ｓｗは、走査重みであり、例えば１６である。ここで、候補にある点の範囲を決めておく。

ここで、Ｎｏであれば、ステップＳ３３−３−８に進む。Ｙｅｓであれば、ステップＳ３３−３−５において、(x,y)と(Xg,Yg)より距離差dを算出する（図１６（ｂ）参照）。ここで、ステップＳ３３−３−６において、d<dmであるか否かを判定する。Ｙｅｓの場合には、ステップＳ３３−３−７に進み、x1=x, y１=y, dm=d（前極側エッジ点候補）と更新して、Ｎｏの場合とともに、ステップＳ３３−３−８に進む。このような処理を繰り返すことにより、前極側のエッジ点が決まる。

ステップＳ３３−３−８では、ｘを１だけインクリメントする。次いで、ステップＳ３３−３−９において、ｘとＸｓとを比較し、x<Xsの場合には、ステップＳ３３−３−３に戻る。x≧Xsの場合には、ステップＳ３３−３−１０に進み、ｙとＹｓとを比較し、y≧Ysの場合には、ステップＳ３３−３−３に戻る。y≧Ysの場合には、ステップＳ３３−３−１２に進みX1←x1, Y1←y1として、前極側エッジ点の補正処理を終了する。

図３２は、図３１のステップＳ３３−３−１およびＳ３３−３−５の処理である前極側エッジ点と重心点方向の距離差算出処理の流れを示すフローチャート図である（図１６、ステップＳ３３−４）。尚、図３２では本処理を前極側エッジ点の補正だけでなく、後述する重心点の補正にも流用するため、前極側エッジ点に対して前方エッジ点、重心点に対して後方エッジ点という汎用化した表記にしている。図３１のステップＳ３３−３−１の(X1,Y1)およびＳ３３−３−５の(x,y)はいずれも、図３２の前方エッジ点(Xs,Ys)に対応し、図３１のステップＳ３３−３−１およびＳ３３−３−５の(Xg,Yg)はいずれも、図３２の後方エッジ点(Xe,Ye)に対応する。距離差の算出は、まず、ステップＳ３３−４’−１において、
前方エッジ点(Xs,Ys)，後方エッジ点(Xe,Ye)とし、
前方エッジ点から後方エッジ点の方向の単位ベクトル
dx=(Xe-Xs)/L, dy=(Ye-Ys)/Lを求める。

但し、L={(Xe-Xs)²+ (Ye-Ys)²}^1/2
前極側エッジ点Ｐ_ｆからの距離ｒ（初期値：2, 上限：12）を定義し、
出力距離差の初期値=0として、ｒ＝２から＝１２まで１１回計算する。

ステップＳ３３−４’−２において、
xc=dx×r+Xs, yc=dy×r+Ysの計測用中間点を規定し（前方エッジ点から後方エッジ点に向かう軸上に規定される計測用中間点 (xc,yc)が求まる。）、前記計測用中間点より前方エッジ点から後方エッジ点に向かう軸と直交する上方向に直線を延ばし、眼球領域の上側の輪郭上の交点(xc1,yc1)を以下のように算出する。

ステップＳ３３−４’−３において、xc1=yc1=-1, t=0.5（０．５画素単位）を代入して変数を初期化し、ステップＳ３３−４’−４において、
x= -dy×t+xc, y=dx×t+ycで表される直線、すなわち、前方エッジ点から後方エッジ点に向かう軸と直交する方向の単位ベクトルを、ｔという値で変化させて上側に移動させる。ｔを増していきながら線を上方向に延ばして輪郭とぶつかったところを探す、すなわち、ステップＳ３３−４’−５に示すように、
x<Xs, y<Ys, ImageM(x,y)=1であるか否かを判定することで見つける。

Ｙｅｓの場合には、ステップＳ３３−４’−６において、上方向に延ばした線と輪郭とが交わった点を（図１６（ｂ））
xc1=x, yc1=yとし、
Ｎｏの場合には、ステップＳ３３−４’−７でｔに０．５を足してステップＳ３３−４’−４に戻る。

次は、前記計測用中間点より前方エッジ点から後方エッジ点に向かう軸と直交する下方向に直線を延ばし、眼球領域の下側の輪郭上の交点(xc2,yc2)を同様に算出する。ステップＳ３３−４’−６からステップＳ３３−４’−８に進み、xc2=yc2=-1, t=-0.5とする。次いで、ステップＳ３３−４’−９において、
x= -dy×t+xc, y=dx×t+ycとする。

次いで、ステップＳ３３−４’−１０において、
x<Xs, y<YsでImageM(x,y)=1であるか否かを判定し、Ｎｏの場合には、ステップＳ３３−４’−１２においてｔから０．５を減算し、ステップＳ３３−４’−９に戻る。Ｙｅｓの場合には、ステップＳ３３−４’−１１において、延ばした線と輪郭とが交わった点を（図１６（ｂ））
xc2=x, yc2=y
として、ステップＳ３３−４’−１３において、前記計測用中間点より上側の輪郭上の交点までの距離と前記計測用中間点より下側の輪郭上の交点までの距離との距離差を以下のように算出して距離差の総和値ｄに加算する。
ｄ←ｄ+|(xc1-xc)²+(yc1-yc)²-(xc2-xc)²-(yc2-yc)²|
とし（ここで、求めたｄは、点(xc,yc)を基準にして、(xc1,yc1)と(xc2,yc2)との間の距離、の総和である。）、ステップＳ３３−４’−１４でｒをセンター方向（重心点）に動かし、ステップＳ３３−４’−１５で
r<上限値(12)
であるか否かを判定することで、計測用中間点(xc,yc)を動かしながら複数個所にわたって上下方向の距離差を求めて、その総和をｄとし、ｄが出力される。

Ｙｅｓの場合には、ステップＳ３３−４’−２に戻り、Ｎｏの場合には、処理を終了する（Ｅｎｄ）。この処理により、前極側エッジ点と重心点方向の距離差ｄを算出することができる。

図３３は、重心点の補正と決定処理（Ｓ３３−４）の流れを示すフローチャート図である（図１６（ｃ）参照）。基本的な考え方は図３１と同様であり、図３１は重心点を仮決めし、仮の重心点から前極側エッジ点を求めた。そして、距離差ｄが最も小さくなるときが前極側エッジ点となる。異なる点は、図３１では輪郭線上を通るのに対して、重心は内部で変動する点が異なる。ステップＳ３３−４−１において、(Xg,Yg)と(X1,Y1)より距離差dを算出する。(Xg,Yg)はいずれも、図３２の前方エッジ点(Xs,Ys)に対応し、(X1,Y1)は図３２の後方エッジ点(Xe,Ye)に対応させ、同様なロジックで距離差ｄを算出する。ここで、ステップＳ３３−４−２において、重心点(xg,yg)、初期値xg=Xg, yg=Yg, 最小距離差dm=dを設定する。次いで、ステップＳ３３−４−３において、左右方向であるか否かを判定する。Ｎｏの場合（上下方向）には、ステップＳ３３−４−１１に進み、Ｙｅｓの場合（左右方向）にはステップＳ３３−４−４に進む。

ステップＳ３３−４−４（左右方向）においては、重心をｙ軸方向に移動させる（ｘ軸方向は動かさない（ｘｇは初期値のまま変更しない）。
Ymin={Yg*(100-Sw)+Y1*Sw}/100
Ymax={Yg*(100-Sw)+Y2*Sw}/100, y=Ymin
とする。これにより、動かす範囲を決めて、不要な範囲を動かさないようにする。ここで、Ｓｗ：走査重みは例えば３０である。

次いで、ステップＳ３３−４−５において、Ｙ方向に関して、(Xg,y)と(X1,Y1)より距離差dを算出する。(Xg,y)はいずれも、図３２の前方エッジ点(Xs,Ys)に対応し、(X1,Y1)は図３２の後方エッジ点(Xe,Ye)に対応させ、同様なロジックで距離差ｄを算出する。

まず、ステップＳ３３−４−６において、ｄ＜dmであるか否かを判定する。Ｙｅｓの場合には、ステップＳ３３−４−７において、
yg=y, dm=d（重心候補）
とし、Ｎｏの場合とともに、ステップＳ３３−４−８において
y←y+1
とし、ステップＳ３３−４−９において、y>Ymaxを判定する。

Ｎｏの場合には、ステップＳ３３−４−６に戻り、Ｙｅｓの場合には、ステップＳ３３−４−１０において
Yg=ygとする。

一方、ステップＳ３３−４−１１において、Ｓｗは走査重み（例３０）である。これにより、上下方向の動かす範囲を決めて、不要な範囲を動かさないようにする。
Xmin={Xg*(100-Sw)+X1*Sw}/100
Xmax={Xg*(100-Sw)+X2*Sw}/100, x=Xmin

ステップＳ３３−４−１２において、
Ｘ方向に関して、(x,Yg)と(X1,Y1)より距離差dを算出する。(x,Yg)はいずれも、図３２の前方エッジ点(Xs,Ys)に対応し、(X1,Y1)は図３２の後方エッジ点(Xe,Ye)に対応させ、同様なロジックで距離差ｄを算出する。

次いで、ステップＳ３３−４−１３において、d<dmを判定する。Ｎｏの場合には、ステップＳ３３−４−１３に戻り、Ｙｅｓの場合には、ステップＳ３３−４−１７に進み、Xg＝xgとすることでＸｇが求まる。
これにより、図１６（ｃ）に示すように、重心点の補正と決定処理が終了する。

次に、図３４は、後極側エッジ点の補正と眼軸の決定処理の流れを示すフローチャート図である（図１７の眼軸の自動設定処理Ｓ３３）。Ｐｆ、Ｐｇ’が決まったため、次いで、Ｐｆ、Ｐｇ’の延長線上に来るようにＰｂを求める。まず、ステップＳ３５−５−１において、(X1,Y1)-(Xg,Yg)-(X2,Y2)より直線性dの判定処理を行う。ステップＳ３５−５−２において、後極側エッジ点(x2,y2)、初期値x2＝X2, y2＝Y2、最小判定値dm＝d、 x＝y＝0を設定する。ステップＳ３３−５−３において、ImageM(x,y)と1とを比べる。ImageM(x,y)が１でない場合には、ステップＳ３３−５−８に進み、ImageM(x,y)が１の場合には、ステップＳ３３−５−４に進む。ここで、ｘ、ｙの値が以下の条件で示される所定の範囲に含まれるか否かを確認する。
左方向：x>{X2*(Sw-1)+Xg}/Sw
右方向：x<{X2*(Sw-1)+Xg}/Sw
上方向：y>{Y2*(Sw-1)+Yg}/Sw
下方向：y<{Y2*(Sw-1)+Yg}/Sw
を満たすか？
を確認する（Ｓｗ： 走査重み（例は４））。

Ｎｏの場合には、ステップＳ３３−５−８に進む。Ｙｅｓの場合には、ステップＳ３３−５−５において、(X1,Y1)-(Xg,Yg)-(x,y)より直線性dを判定する処理を行う。

まず、制約を決めるため、ステップＳ３３−５−６において、d＜dmを判定する。Ｎｏの場合には、ステップＳ３３−５−８に進む。Ｙｅｓの場合には、
x2＝x, y2＝y, dm＝d（後極側エッジ点候補）
となる。次いで、ステップＳ３３−５−８において、
x←x+1
とし、ステップＳ３３−５−９において、xとXsとを比較する。ここで、x<Xsの場合には、ステップＳ３３−５−３に戻る。x≧Xsの場合には、
x=0, y←y+1として、ステップＳ３３−５−１１においてyとYｓとを比較する。y≧Ysの場合には、ステップＳ３３−５−３に戻り、y≧Ysの場合にはステップＳ３３−５−１２において、
X2←x2, Y2←y2
として、処理を終了する（Ｅｎｄ）。これにより、Ｘ２、Ｙ２を補正することができる（図１７（ｂ）参照）。

図３５は、図３４のＳ３３-５-１と、Ｓ３３-５-５の処理内容の具体的な例である３頂点の直線性の判定処理の流れを示すフローチャート図である（図１７（ｂ）参照）。図３４のＳ３３-５-１における(X1,Y1)は図３５における(xs,ys)に、(Xg,Yg)は図３５における(xm,ym)に、(X2,Y2)は図３５における(xe,ye)に対応し、Ｓ３３-５-５における(X1,Y1)は図３５における(xs,ys)に、(Xg,Yg)は図３５における(xm,ym)に、(x,y)は図３５における(xe,ye)に対応する。まず、ステップＳ３３−５−１−１において、与えられた３頂点を(xs,ys), (xm,ym), (xe,ye)とする。また、dx=xe-xs, dy=ye-ys、出力判定値d=0とする。

ステップＳ３３−５−１−２において、|dx|>|dy|を判定する。Ｎｏの場合（上下方向）には、ステップＳ３３−５−１−４に、Ｙｅｓの場合（左右方向）にはステップＳ３３−５−３に進む。

ステップＳ３３−５−１−３において、ｄ=|dy×(xm-1)/dx+ys-ym|により、ｄを求め、ステップＳ３３−５−１−４において、ｄ=|dx×(ym-1)/dy+xs-xm|によりｄを求め、処理を終了する（Ｅｎｄ：図１７（ｂ）参照））。
ここまでが、眼軸の設定処理である。

以降は、手動の調整を行う。
図３６は、手動設定エッジ点の輪郭線上への自動補正処理（Ｓ３５−３）を示すフローチャート図である。図１８の輪郭線上の点を見つける処理である。

まず、ステップＳ３５−３−１において、
補正対象の前極側エッジ点を(x1,y1)、補正後の前極側エッジ点を(X1,Y1)、重心点を(Xg,Yg)とし、後極側エッジ点(X2,Y2)は固定とする。尚、具体的な説明は省略するが、逆に、前極側エッジ点(X1,Y1)を固定とし、補正対象の後極側エッジ点を(x2,y2)、補正後の後極側エッジ点を(X2,Y2)、重心点を(Xg,Yg)とする方法も同様に実現できる。補正前の前極側エッジ点を用いて眼軸方向の単位ベクトルを以下のように定義する。
L={(x1-X2)²+ (y1-Y2)²}^1/2
dx=(x1-X2)/L, dy=(y1-Y2)/L
t=0.5
補正後の前極側エッジ点を(x,y)とし、初期値をｘ＝x１、ｙ＝y１とする。

ステップＳ３５−３−２において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、ImageM(ｘ,ｙ)=1であるか否かを判定する。Ｙｅｓの場合には、ステップＳ３５−３−１０に進み、Ｎｏの場合には、ステップＳ３５−３−２に進み、
x=dx×t+x1, y=dy×t+y1
と変更する。次いで、ステップＳ３５−３−４において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、
ImageM(ｘ,ｙ)=1
であるか否かを判定し、Ｙｅｓの場合にはステップＳ３５−３−１０に進み、Ｎｏの場合には、ステップＳ３５−３−に進み、
x=-dx×t+x1, y=-dy×t+y1
と変更する。次いで、ステップＳ３５−３−７において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、ImageM(ｘ,ｙ)=1であるか否かを判定する。Ｙｅｓの場合には、ステップＳ３５−３−１０に進み、Ｎｏの場合には、ステップＳ３５−３−８に進み、t←t+0.5とする。次いで、ステップＳ３５−３−９において、
t<L/2であるか否かを判定し、
Ｎｏの場合には、ステップＳ３５−３−３に戻り、Ｙｅｓの場合には、Ｅｒｒｏｒ（補正不可）となる。ステップＳ３５−３−１０においては、
X1=x, Y1=y
Xg=(X1+X2)/2
Yg=(Y1+Y2)/2
として、処理を終了する（Ｅｎｄ）。これにより、手動設定されたエッジ点の輪郭線上への自動補正処理が行われる。

図３７は、既に決定された眼軸と直交する方向に９０度および−９０度エッジ点の追加を行う処理（ステップＳ３８−１、図１９（ａ）参照）の流れを示すフローチャート図である。まず、ステップＳ３８−１−１において、前極側エッジ点を基点にした後極方向への眼軸上の単位ベクトルを求める。図３６までのＸｇ，Ｙｇは、眼軸の自動設定を行うために導入した眼球の幾何学的な重心であり、図３７、図３８のＸｏ，Ｙｏは、計測用重心であり、人為的に設定した眼軸上の後極点から一定距離離れた値であるため、両者では値がかわってくる（ただし、正常眼ではほぼ一致するように計測用重心を設定する）。眼軸上の単位ベクトルは、
dx=(X2-X1)/L, dy=(Y2-Y1)/L
ただし、L={(X2-X1)²+ (Y2-Y1)²}^1/2
とする。

ステップＳ３８−１−２において、ｔ＝０．５とする。ステップＳ３８−１−３において、計測用重心より上記単位ベクトルの方向に伸ばした点として、
x=-dy×t+X0, y=dx×t+Y0
を定義する。

次いで、ステップＳ３８−１−４において、上方向にスキャンし、輪郭線にぶつかる点を求める。対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、
x<Xs, y<Ys, ImageM(x,y)=1
を判定し、Ｎｏの場合には、ステップＳ３８−１−５において、
t¬t+0.5
とし、ステップＳ３８−１−３に戻る。

Ｙｅｓの場合には、ステップＳ３８−１−６に進み、
Xp9=x, Yp9=y
とする。
ステップＳ３８−１−７において、ｔ＝０．５とする。

次いで、ステップＳ３８−１−８において、下方向にスキャンし、輪郭線にぶつかる点を求める。
x= -dy×t+X0, y=dx×t+Y0とし、
ステップＳ３８−１−９において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、
x<Xs, y<Ys, ImageM(x,y)=1であるか否かを判定する。Ｎｏの場合には、ステップＳ３８−１−１０に進み、t¬t-0.5とした後、ステップＳ３８−１−８に戻る。

Ｙｅｓの場合には、ステップＳ３８−１−１１において、
Xm9=x, Ym9=y
とする。
これにより計測用エッジ点を追加することができる。

図３８は、θ_ｗ度および−θ_ｗ度の計測用エッジ点（±第１の角度）の追加処理の流れを示すフローチャート図である（ステップＳ３８−２、図１９（ａ）、（ｃ）参照）。ここでは、θ_ｗが４５度の例を示すが、第１の角度は、９０度より狭く、後述する第２の角度よりも広い角度である（４０度から９０度）。

まず、ステップＳ３８−２−１（プラス４５度の計測用エッジ点）で、θ_ｗ度方向の単位ベクトルを計算する。本実施例では、以下ｄｘ、ｄｙは４５度方向の単位ベクトルである。
dx={(X2-X0)/L+(Xp9-X0)/Lp}/2
dy={(Y2-Y0)/L+(Yp9-Y0)/Lp}/2
ただし、L={(X2-X0)²+ (Y2-Y0)²}^1/2、Lp={(Xp9-X0)²+ (Yp9-Y0)²}^1/2である。

次いで、ステップＳ３８−２−２で、ｔ＝０．５とする。ステップＳ３８−２−３で、計測用重心より上記単位ベクトルの方向に伸ばした点として
x=dx×t+X0, y=dy×t+Y0を定義する。次いで、ステップＳ３８−２−４において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、x<Xs, y<Ys, ImageM(x,y)=1であるか否かを判定する。ここでＮｏであれば、ステップＳ３８−２−５において、ｔに０．５を加算する。Ｙｅｓの場合には、ステップＳ３８−２−６において、Xp4=x, Yp4=yと計測用エッジ点が求める（Ｅｎｄ）。

また、ステップＳ３８−２−１１において（マイナスθ_ｗの計測用エッジ点）、−θ_ｗ度方向の単位ベクトルを計算する。本実施例では、以下ｄｘ、ｄｙは−４５度方向の単位ベクトルである。
dx={(X2-X0)/L+(Xm9-X0)/Lm}/2
dy={(Y2-Y0)/L+(Ym9-Y0)/Lm}/2
ただし、L={(X2-X0)²+ (Y2-Y0)²}^1/2 、Lm={(Xm9-X0)²+ (Ym9-Y0)²}^1/2である。
次いで、ステップＳ３８−２−１２で、ｔ＝０．５とする。ステップＳ３８−２−１３で、計測用重心より上記単位ベクトルの方向に伸ばした点として
x=dx×t+X0, y=dy×t+Y0を定義する。次いで、ステップＳ３８−２−１４において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、x<Xs, y<Ys, ImageM(x,y)=1であるか否かを判定する。ここでＮｏであれば、ステップＳ３８−２−１５において、ｔに０．５を加算する。Ｙｅｓの場合には、ステップＳ３８−２−１６において、Xｍp4=x, Yｍ4=yと計測用エッジ点が求める（Ｅｎｄ）。
以上により、θ_ｗ度及び−θ_ｗ度の計測用エッジ点を追加することができる。

図３９は、θ_ｎ度および−θ_ｎ度の計測用エッジ点（±第２の角度）、の追加処理の流れを示すフローチャート図である（ステップＳ３８−３、図１９（ａ）、（ｄ）参照）。ここでは、第２の角度が２２．５度の例を示すが、第２の角度は、上記第１の角度よりも狭い角度である（１０度から４０度）。

ステップＳ３８−３−１において、以下のように、＋２２．５度方向の単位ベクトルを定義する。
dx={(X2-X0)/L+(Xp4-X0)/Lp}/2
dy={(Y2-Y0)/L+(Yp4-Y0)/Lp}/2
ただし、L={(X2-X0)²+ (Y2-Y0)²}^1/2 、Lm={(Xm4-X0)²+ (Ym4-Y0)²}^1/2である。

次いで、ステップＳ３８−３−２で、ｔ＝０．５とする。ステップＳ３８−３−３で、計測用重心より上記単位ベクトルの方向に伸ばした点として
x=dx×t+X0, y=dy×t+Y0を定義する。次いで、ステップＳ３８−３−４において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、x<Xs, y<Ys, ImageM(x,y)=1であるか否かを判定する。ここでＮｏであれば、ステップＳ３８−３−５において、ｔに０．５を加算する。Ｙｅｓの場合には、ステップＳ３８−３−６において、Xp２=x, Yp２=yとエッジ点が求める（Ｅｎｄ）。

また、ステップＳ３８−３−１１において、以下のように、−２２．５度方向の単位ベクトルを定義する。
dx={(X2-X0)/L+(Xm4-X0)/Lm}/2
dy={(Y2-Y0)/L+(Ym4-Y0)/Lm}/2
ただし、L={(X2-X0)²+ (Y2-Y0)²}^1/2 、Lm={(Xm4-X0)²+ (Ym4-Y0)²}^1/2である。

次いで、ステップＳ３８−３−１２で、ｔ＝０．５とする。ステップＳ３８−３−１３で、計測用重心より上記単位ベクトルの方向に伸ばした点として
x=dx×t+X0, y=dy×t+Y0を定義する。次いで、ステップＳ３８−３−１４において、対象画素（ｘ、ｙ）が輪郭上にあるか否か、即ち、x<Xs, y<Ys, ImageM(x,y)=1であるか否かを判定する。ここでＮｏであれば、ステップＳ３８−３−１５において、ｔに０．５を加算する。Ｙｅｓの場合には、ステップＳ３８−３−１６において、Xｍ２=x, Yｍ２=yと計測用エッジ点が求める（Ｅｎｄ）。

次に、後極部における円弧面積の算出方法を説明する。面積比率は、−θ_ｗ度から＋θ_ｗ度の範囲の広領域での算出と、−θ_ｎ度から＋θ_ｎ度の範囲の狭領域での算出の２段階で行われるが、始めに、広領域での算出方法について説明する。

図４０（ａ）は左方向投影画像の場合の図であり、図４０（ｂ）は、求めた広領域での面積値を基に比率を算出する方法を示す表である。

図４１は、図４０で説明した後極部における円弧面積の算出処理例を示すフローチャート図である。まず、ステップＳ４０ｃ−１において、上部面積Ap＝０、ｖ＝０と初期化し、ステップＳ４０ｃ−２でｕ＝０と初期化する。次いで、ステップＳ４０ｃ−３において、
x=(Xp4-X0)×u/Lup+(X2-X0)×v/Lv
y=(Yp4-Y0)×u/Lup+(Y2-Y0)×v/Lv
の式により、＋４５度方向と０度方向に囲まれた領域内のｘとｙを求める。

ステップＳ４０ｃ−４において、（ｘ、ｙ）が面積の算出対象領域内にあるか否か、即ち、ImageM(x,y):0を判定する。

０でない場合には、同一画素が既に重複してカウントされていないかをBuff(x,y):0で確認し、Buff(x,y)＝0の場合のみBuff(x,y)=1、Ap¬Ap+1を代入して、ステップＳ４０ｃ−７に進み、u¬u+0.5とする。ステップＳ４０ｃ−８で、uとLupの大小を比較し、u<Lupであれば、ステップＳ４０ｃ−３に戻り、u³Lupであれば、ステップＳ４０ｃ−９に進む。ｖについて判定をし、その結果、u³Lupであれば、上部面積Ａｐが決定される。同様に処理によりＳ４０ｃ−１１以降の処理でＡｍについて求める。Ａｍの算出において、Ａｐの算出と異なる点は、Ｓ４０ｃ−１５、Ｓ４０ｃ−１６の処理で、Ａｐの算出において同一画素が既に重複してカウントされている場合は、ＡｐとＡｍの双方に重複する境界領域であることを示すため、Ａｍ側で二重カウントを抑止するだけでなく、Ｓ４０ｃ−１６の処理で、Ａｐ側に既にカウントされた当該１画素を取り消す処理を行う。このようにして、図２０に示す後極部における円弧面積を求めることができる。

続いて狭領域での算出方法について説明する。図４２（ａ）は、後極部における円弧面積の算出方法（狭領域）を示す図である。この図は、左方向投影画像の場合の例である。上記で求めたＸｐ２、Ｙｐ２と、Ｘｍ２、Ｙｍ２を用いて算出を行う。図４２（ｂ）に示すように、左方向画像、右方向画像、上方向画像、下方向画像におけるＥ／Ｎ比率、Ｔ／Ｂ比率を求めることができる。

図４３は、後極部における円弧面積の算出方法（狭領域）を示すフローチャート図である。前述した広領域と同様に計算できる。単位ベクトルは、上記において定義したものである（数１）。

ここで、＋２２．５度方向、−２２．５度方向、０度方向に対して、単位ベクトルUp→=(Xp2-X0,Yp2-0Y)/Lup、Um→=(Xm2-X0,Ym2-0Y)/Lum、V→=(X2-X0,Y2-Y0)/Lvを求める。
ただし、Lup={(Xp2-X0)²+ (Yp2-Y0)²}^1/2、Lum={(Xm2-X0)²+ (Ym2-Y0)²}^1/2、Lv={(X2-X0)²+ (Y2-Y0)²}^1/2

輪郭抽出後の投影画像データをImageM(x,y)={0:背景，1:輪郭，2:内部} (x=0,…,Xs; y=0,…,Ys) Xs=Ys=300 [pixel]とし、併せて同サイズの画像バッファBuff(x,y)を定義し、画像バッファの初期値は全て０にする。画像バッファBuff(x,y)は、同じ画素を２回カウントするのを防止する。

上部面積をAp2、下部面積をAm2とし、初期値は各々０とする。算出されたAp2とAm2を用いて、面積比率を算出する。

まず、ステップＳ３９−１において、Ap2=0, v=0を設定し、ステップＳ３９−２において、ｕ＝０とする。ｖ，ｕは媒介変数である。ステップＳ３９−３において図２０（ａ）の後極部における左右の上側の円弧面積の算出を以下の式により計算する。＋２２．５度方向と０度方向に囲まれた領域内のｘとｙを求める。
x= (Xp2-X0)×u/Lup+ (X2-X0)×v/Lv
y= (Yp2-Y0)×u/Lup+ (Y2-Y0)×v/Lv

ステップＳ３９−４で、（ｘ、ｙ）が面積の算出対象領域内にあるか否か、即ち、ImageM(x,y):0であるか否かを計算する。Ｎｏの場合には、ステップＳ３９−５で同一画素が既に重複してカウントされていないかをBuff(x,y):0で確認し、０であれば、ステップＳ３９−６でバッファに同じ値Buff(x,y)=1を格納し、左右方向上側の面積Ap2に１を加算する。そして、ステップＳ３９−４でImageM(x,y)＝0である場合とともに、ステップＳ３９−８において、ｕに０．５を加算し、ステップＳ３９−８でｕとＬｕｐとを比較する。ここで、ｕ＜Ｌｕｐであれば、ステップＳ３９−２に戻り、ｕ³Ｌｕｐであれば、ステップＳ３９−９に進み、ｖに０．５を加算する。そして、ステップＳ３９−１０に進み、ｖとＬｖとを比較し、ｖ＜ＬｖであればステップＳ３９−２に戻り、ｖ³Ｌｖであれば、ステップＳ３９−１１において、Ａｍ２＝０、ｖ＝０とし、ステップＳ３９−１２でｕ＝０とする。これで、Ａｐ２が求まる。

次いで、ステップＳ３９−１３において、−２２．５度方向と０度方向に囲まれた領域内のｘとｙを求める。
x= (Xm2-X0)×u/Lum+ (X2-X0) ×v/Lv
y= (Ym2-Y0)×u/Lum+ (Y2-Y0)×v/Lv
により、図２０（ａ）の後極部における左右の下側の円弧面積の算出を行う。

次いで、ステップＳ３９−１４において、（ｘ、ｙ）が面積の算出対象領域内にあるか否か、即ち、ImageM(x,y)が0か否かを判断し、ImageM(x,y)＞０であれば、ステップＳ３９−１５に進み、同一画素が既に重複してカウントされていないかをBuff(x,y):0で確認する。Buff(x,y)が0であれば、ステップＳ３９−１７に進み、
Buff(x,y)=2
Am2←Am2+1
とし、
Buff(x,y)が0でなければ、ステップＳ３９−１６に進み、
Buff(x,y)=2
Ap2←Ap2-1
とする。これで、上側において既にカウントされる境界部分のカウントを除く。

そして、ステップＳ３９−１８でｕに０．５加算し、ステップＳ３９−１９でｕとＬｕｍとを比較する。ｕ＜Ｌｕｍであれば、ステップＳ３９−１２に戻り、ｕ³Ｌｕｍであれば、ステップＳ３９−２０でｖに０．５を加算する。そして、ステップＳ３９−２１において、ｖとＬｖとを比較し、ｖ＜ＬｖであればステップＳ３９−１２に戻る。ｖ³Ｌｖであれば、処理を終了する。これで、Ａｍ２が求まる。これにより、図２０の後極部における円弧面積をＡｐ２とＡｍ２として重複しないように画素をカウントして算出することができる。

図４４は、後極部３エッジ点（尖鋭度）の成す角の算出方法の一例を示す図である。

尚、前述の面積比率の算出と異なり、上記なす角の算出は、狭領域の面積に対応する領域のみで行う。
以上で、解析処理のアルゴリズムについての説明を行った。

次に、実際に得られる投影画像及び測定結果について説明する。
図４５は、解析対象の６方向投影画像例であり、左右、上下、前後の６方向からの投影画像が得られていることがわかる。図４６は、解析における解析パラメータの設定画面の例を示す図である。設定パラメータとしては、二値化しきい値（０〜２５５の範囲、通常の設定値：２５）、削除スポット面積（単位：画素、通常の設定値：２００）、探索範囲重み（前極、後極、通常の設定値：１６と４）、水晶体エッジの測定範囲（図１５Ａ（２）（３）においてＰｃを移動させる範囲、通常の設定値：２〜１２）、重心探索範囲（図１５Ａ（３）においてＰｇを移動させる範囲、単位％、通常の設定値：３０）、計測の原点（＝計測用重心、後極側エッジ点からの距離、通常の設定値：８７画素）、対称性算出角度（広領域の面積算出角度範囲、通常の設定値：４５度）、尖鋭度算出角度（尖鋭度算出および狭領域の面積算出角度範囲、通常の設定値：２２度）、解析眼球が左眼であるか右眼であるか、標準眼軸長上限（単位：画素、通常の設定値：１７５画素）、左右対称性下限値（％、通常の設定値：９０％）、上限値（％、通常の設定値：１１０％）、上下対称性下限値（％、通常の設定値：９０％）、上限値（％、通常の設定値：１１０％）、後極尖鋭度下限値（通常の設定値：１４０度）、上限値（通常の設定値：１５０度）を設定することができるが、一般的には、解析眼球が左眼であるか右眼であるのみを設定すれば良い。

図４７は、解析結果の一例を示す図である。上記解析処理に基づく解析結果を一画面で見ることができる。例えば、Ｌ（左方向投影像）、Ｒ（右方向投影像）、Ｔ（上方向投影像）、Ｂ（下方向投影像）について、眼軸長（単位：画素）、後極尖鋭度（０〜３６０度）、狭領域面積比率（％）、広領域面積比率（％）、面積比率（Ｔ／Ｂ、Ｅ／Ｎ、狭領域または広領域面積比率のいずれかの値、１００との差が大きい値）、後極ピーク位置（座標値、画素単位）、左右方向の対称性（耳寄り、左右対称、鼻寄りのいずれかの判定結果）、上下方向の対称性（上寄り、左右対称、下寄りのいずれかの判定結果）、後極尖鋭タイプ（樽型、中間、紡錘型のいずれかの判定結果）、タイプ分類（表１の４桁）が表示されるようになっている。タイプ分類は２０００であり、後極尖鋭タイプが２の樽型であり、その他の眼軸長、左右方向の対称性、上下方向の対称性は正常である。この画面より、解析結果と、病態とが一目でわかるようになっている。

図４８は、図４５に対応する図であり、解析時の補助線設定状況の確認画面の一例を示す図である。図４８においては、上下左右前後の眼球の投影画像において、眼軸について修正を行うための補助線（計測用重心から各種計測用エッジ点を結ぶ直線）を示している。この確認画面において、眼軸を修正し、指示した補助線で再解析を行うことができる。

以上は、正常な場合の例であるが、以下に、強度近視の場合の例について示す。
図４９は、強度近視の場合の解析対象の６方向投影画像の例を示す図である。図４５に示す画像とは、形状がかなり異なっていることがわかる。図５０は、解析パラメータの設定画面例を示す図であり、この画面は、図４６と同様である。同じ設定において、解析を行った結果の例について図５１を例にして説明する。図４６とは、解析結果がかなり異なっていることがわかる。特に、眼軸長、後極尖鋭度、左右方向の対称性、上下方向の対称性、後極尖鋭度のタイプなどが異なっている。タイプ分類において、表１において、タイプ分類が１１であり、眼軸長が異常、左右方向の対称性が鼻寄り（９０％未満）であることがわかる。図５２は、解析時の補助線設定状況の確認画面が示される。図５２においては、上下左右前後の眼球の投影画像において、眼軸について修正を行うための補助線を示している。この確認画面において、眼軸を修正し、指示した補助線で再解析を行うことができる。

図５３は、解析対象の６方向投影画像（軽症）の６方向投影画像の例を示す図である。図４９と比較すると、やや図４５に近い投影画像になっていることがわかる。この場合には、図５４に示すように、左眼を解析対象としている。その結果が、図５５であり、タイプ分類をみると、１２１であり、眼軸長が異常、左右方向の対称性が耳寄り（１１０％以上）、上下方向の対称性が下寄り（９０％未満）であることがわかる。

図５６は、解析時の補助線設定状況の確認画面の例を示す図である。図５６に示すように、眼軸自動設定については、画面の平行・垂直方向とはずれており、明らかに不適切であることがわかる。図５７は、解析時の補助線設定状況の確認画面であって、眼軸の修正（この例では手動修正）を付加した様子を示す図である。図５７に示すように、眼軸の補助線は、画面の平行・垂直方向と一致しており、正しく修正されていることが確認できる。

図５８は、解析パタメータの設定画面であり、例えば図５４に示す設定画面から設定値が変更されていないが、設定値を画面で確認することができる。

図５９は、修正後の再解析結果の一例を示す図である。図５５と比較して、左右方向の対称性が耳寄りから左右対称に変化しており、後極尖鋭度は中間型から紡錘型に変化しており（正常眼は図４７に示される通り樽型）、タイプ分類は１２１から１１０１に変化している。このように、眼軸を変更した後に解析結果を得ることで、正しい値を得ることができる。

上記の実施の形態において、添付図面に図示されている構成等については、これらに限定されるものではなく、本発明の効果を発揮する範囲内で適宜変更することが可能である。その他、本発明の目的の範囲を逸脱しない限りにおいて適宜変更して実施することが可能である。上記眼球部ＭＲＩの画像処理方法をコンピュータに実行させるためのプログラム、当該プログラムを記録するコンピュータ読み取り可能な記録媒体であっても良い。

（まとめ）
本発明によれば、患者に負担をかけずに「強度近視」の早期診断に寄与でき、読影医の主観に依存しない客観性のある診断データを提供することができる。

また、安価な汎用パーソナルコンピュータのソフトウェアで実現することができ、ＧＰＵ（Graphics Processing Unit）など特殊なグラフィックスボードの増設も必要とないため、安価な診断サービスを提供できる。

また、画像入力系も新規に眼科専用設備を導入することなく、既存のポピュラーなＭＲＩ機器を流用し、放射線を使用しないので被爆の問題が無く、定期的に繰り返し診断を行うことができ、眼科向けのカスタマイズも必要としない。即ち、本願発明は、新規なハードウェア設備の投資をすることなく、安価のシステムを構築することができる。

（付記）
本発明は、以下の発明を含む。
（１）眼科診断機器による診断を行う第１のステップと、
前記第１のステップで診断ができない場合に、頭部のＭＲＩ撮影及び眼球部位の３次元ボクセル画像を抽出するステップと、
前記眼球部位の３次元ボクセル画像を基に、上下／左右／前後の６方向のいずれかの方向の投影画像データを左右眼球別にレンダリングするステップと、
前記投影画像に対して、画像解析による眼球の形状パラメータである、眼軸長、後極先鋭度、左右対称性、上下対称性、後極突出数の少なくともいずれか１のパラメータに基づいて、強度近視を診断するステップと、
を有することを特徴とする強度近視の診断方法。
（２）前記形状パラメータ毎に、強度近視のタイプの分類処理を行うステップを有することを特徴とする（１）に記載の強度近視の診断方法。

このように、一般的な眼科診断で強度近視の診断が確定しない場合に、本実施の形態による、頭部ＭＲＩ撮影データを用いた眼科疾患の画像解析技術を用いると良い。

本発明は、眼科疾患の画像解析装置に利用できる。

Ａ…画像解析装置、１…画像入力部、３…ボクセル画像抽出部、５…投影画像算出部、７…投影画像解析部、１１…出力部、７−１…非対称性算出部、７−２…先鋭度合い算出部、７−３…後極突出部ピーク算出部、７−４…眼軸自動設定部、７−５…頂点補正部。

Claims (17)

1. 眼球を含む頭部の３次元ボクセル画像を取得する画像入力部と、
   入力した前記３次元ボクセル画像より、左右眼球別に眼球領域の少なくともいずれか一
   方の眼の３次元眼球ボクセル画像を抽出するボクセル画像抽出部と、
   抽出した３次元眼球ボクセル画像を、少なくとも上下いずれかの方向に投影した第１の
   投影画像を作成する投影画像算出部と、
   前記第１の投影画像に対して眼球前極部頂点と眼球後極部頂点とを結ぶ眼軸を設定し、少なくとも眼球前極部頂点と眼球後極部頂点との間の眼軸長と、眼軸を境界に後極側を左右に分割した領域をもとに左右方向の非対称性の度合いと、を画像解析のパラメータとして算出する投影画像解析部と、
   を備えることを特徴とする眼科疾患の画像解析装置。
2. 前記３次元ボクセル画像は、ＭＲＩを用いて取得したものであることを特徴とする請求項１に記載の眼科疾患の画像解析装置。
3. 前記投影画像算出部は、
   前記３次元眼球ボクセル画像を、少なくとも左右いずれかの方向に投影した第２の投影画像を作成し、
   前記画像解析部は、前記第２の投影画像に対して眼軸を設定し、少なくとも眼球前極部頂点と眼球後極部頂点との間の眼軸長と眼軸を境界に後極側を上下に分割した領域をもとに上下方向の非対称性の度合いを、画像解析のパラメータとして算出することを特徴とする請求項１又は２に記載の眼科疾患の画像解析装置。
4. 前記画像解析部は、
   前記第１又は第２の投影画像のうちの少なくともいずれかに対して、後極の円弧部曲率の尖鋭度合いを、画像解析のパラメータとして算出することを特徴とする請求項１から３までのいずれか１項に記載の眼科疾患の画像解析装置。
5. 前記画像解析部は、
   前記第１又は第２の少なくともいずれか一方の投影画像について、左右又は上下の少なくともいずれか一方の組の前記パラメータを平均した値を算出結果とすることを特徴とする請求項１から４までのいずれか１項に記載の眼科疾患の画像解析装置。
6. 前記投影画像算出部は、
   前記３次元眼球ボクセル画像に対して、後方向に投影した第３の投影画像を作成し、
   前記画像解析部は、前記第３の投影画像に対して、後極突出部に対応するピーク点の個数を、画像解析のパラメータとして算出することを特徴とする請求項１から５までのいずれか１項に記載の眼科疾患の画像解析装置。
7. 前記画像解析部は、前記第３の投影画像に対して、眼球部の輪郭線抽出を行い、更に眼球部の重心点を設定し、前記算出された各ピーク点の位置を、画像解析のパラメータとして、前記重心点を基準とした極座標である重心からの距離と重心からピーク点方向の角度で算出することを特徴とする請求項６に記載の眼科疾患の画像解析装置。
8. 前記画像解析部は、
   前記投影画像に対して眼球部の輪郭線抽出を行い、前記眼球前極部頂点および前記眼球後極部頂点は前記眼球部の輪郭線上に位置し、かつ前記眼球前極部頂点と前記眼球後極部頂点とを結ぶ直線が前記直線の幾つかの頂点より鉛直方向に左右または上下方向に延長した前記輪郭線との交点までの距離が均等になるように、前記眼球前極部頂点および前記眼球後極部頂点を設定することにより、前記眼軸を自動設定することを特徴とする請求項１から３までのいずれか１項に記載の眼科疾患の画像解析装置。
9. 前記画像解析部は、
   前記投影画像に対して眼球部の輪郭線抽出を行い、更に眼球部の重心点を設定し、前記眼球前極部頂点および前記眼球後極部頂点は前記眼球部の輪郭線上に位置し、かつ前記眼球前極部頂点と前記重心点と前記眼球後極部頂点との３点を結ぶ折れ線ができるだけ直線になるように、前記眼球前極部頂点、前記眼球後極部頂点および前記重心点を設定することにより、前記眼軸を自動設定するようにしていることを特徴とする請求項１から３までのいずれか１項に記載の眼科疾患の画像解析装置。
10. 前記画像解析部は、
    前記第１又は第２の投影画像に対して眼球部の輪郭線抽出を行い、前記眼球前極部頂点または前記眼球後極部頂点のいずれかに対して、画面との対話形式により設定する対話型頂点入力部により設定された前記眼球前極部頂点または前記眼球後極部頂点に対して、前記眼球前極部頂点と前記眼球後極部頂点とを結ぶ直線の延長線上で前記輪郭線との交点に自動的に補正する頂点補正部を備えていることを特徴とする請求項１から３までのいずれか１項に記載の眼科疾患の画像解析装置。
11. 前記画像解析部は、
    前記第１又は第２の投影画像に対して眼球部の輪郭線抽出を行い、更に眼球部に計測用重心点を設定し、計測用重心点より前記眼球後極部頂点を結ぶ基準線に対して、計測用重心点を中心に＋第１の角度および−第１の角度の方向の直線Ａと直線Ｂを設定し、前記直線Ａと前記基準線に囲まれた前記輪郭線の内部の面積Ａと、前記直線Ｂと前記基準線に囲まれた前記輪郭線の内部の面積Ｂを算出し、面積Ａと面積Ｂとの比率を前記非対称性の度合いとして与えることを特徴とする請求項１から３までのいずれか１項に記載の眼科疾患の画像解析装置。
12. 前記画像解析部は、
    前記第１又は第２の投影画像に対して眼球部の輪郭線抽出を行い、更に眼球部に計測用重心点を設定し、計測用重心点より前記眼球後極部頂点を結ぶ基準線に対して、計測用重心点を中心に＋第２の角度および−第２の角度の方向の直線Ａと直線Ｂを設定し、前記直線Ａと前記輪郭線との交点Ａと前記眼球後極部頂点と前記直線Ｂと前記輪郭線との交点Ｂとの３頂点のなす角を前記尖鋭度合いとして与えることを特徴とする請求項４に記載の眼科疾患の画像解析装置。
13. 前記画像解析部は、計測用重心点より前記眼球後極部頂点を結ぶ基準線に対して、計測用重心点を中心に＋第２の角度および−第２の角度の方向の直線Ｃと直線Ｄを追加設定し、前記直線Ｃと前記基準線に囲まれた前記輪郭線の内部の面積Ｃと、前記直線Ｄと前記基準線に囲まれた前記輪郭線の内部の面積Ｄを追加算出し、面積Ｃと面積Ｄとの比率と前記面積Ａと面積Ｂとの比率のいずれか大きい値を前記非対称性の度合いとすることを特徴とする請求項１１に記載の眼科疾患の画像解析装置。
14. 前記画像解析部は、
    前記計測用重心点として前記眼球前極部頂点と前記眼球後極部頂点とを結ぶ眼軸上に、前記眼球後極部頂点より所定の距離だけ離した点を与えることを特徴とする請求項１１から１３までのいずれか１項に記載の眼科疾患の画像解析装置。
15. 前記画像解析部は、
    前記所定の距離として正常者の眼球の平均的な眼軸長の１／２を設定することを特徴とする請求項１４に記載の眼科疾患の画像解析装置。
16. 眼球を含む頭部の３次元ボクセル画像を取得する画像入力ステップと、
    入力した前記３次元ボクセル画像より、左右眼球別に眼球領域の少なくともいずれか一方の眼の３次元眼球ボクセル画像を抽出するボクセル画像抽出ステップと、
    抽出した３次元眼球ボクセル画像を、少なくとも上下いずれかの方向に投影した第１の投影画像を作成する投影画像算出ステップと、
    前記第１の投影画像に対して眼球前極部頂点と眼球後極部頂点とを結ぶ眼軸を設定し、少なくとも眼球前極部頂点と眼球後極部頂点との間の眼軸長と、眼軸を境界に後極側を左右に分割した領域をもとに左右方向の非対称性の度合いと、を画像解析のパラメータとして算出する投影画像解析ステップと、
    を備えることを特徴とする眼科用画像解析方法。
17. コンピュータに、請求項１６に記載の眼科用画像解析方法を実行させるためのプログラム。

Images