楕円に内接する平行四辺形

[そら]

いつもお世話になっております。

私は、現在、楕円の中心座標を求めるプログラムを作成しております。

プログラムの内容と致しましては、楕円に内接する平行四辺形を求める→その平行四辺形の各辺の中心点同士を結んだ線が交わった点（楕円の中心座標）を求めるです。

今回は、この手順を踏んで、楕円の中心座標を求めることが目標となっております。

色々と試みたのですが、目的を達成することができず、先が見えなくなってしまったので、今回投稿させて頂きました。

何かアイディア、アドバイス等頂けると幸いです。

何卒、よろしくお願い致します。

なお、プログラムのイメージを添付致しました。

[M.Shiraishi · Accepted solution]

楕円に内接する平行四辺形でないといけませんか？

まず、楕円は、どのように表されているのでしょうか。

2次元マッピング上で、ピクセル的に楕円のエリア（内部or輪郭）が1で、それ以外が0とかの状態でしょうか。

そうであれば、楕円を現すXY座標の最大最小をとってその半分を得れば、つまり楕円に外接する

長方形の中心を求めれば、幾何的にはそれが楕円の中心です。

楕円を表す点のXY座標の最大と最小を見つければよいだけなので、極めて簡便です。

ただし、0.5pixel単位が最小なのと、1pixelはみ出たような点に引っ張られる欠点があります。

なるべく楕円のエリアor輪郭の全てを使うには、例えば、Yの中心を求めるのに、楕円を

Y座標をかえながらX軸に平行な線でスライスしたときの「楕円の幅」（楕円のエリアのpixel数

又は左右の輪郭の座標の差）の2乗をY座標でプロットすると、2次式になります。

（楕円の式を式変形すれば分かる）

なので、これを2次式でフィットすると、その頂点は、楕円のY中心を表します。

同じことをX中心に対しても行えば、楕円の中心が、楕円の全ての情報を使って

フィッティングで割り出せます。かつ、小数で中心が出てきます。

上記2つを、模式図とviで示してみました。（添付）

また、もし、楕円の中が全て1でその外が全て0、とかいう風にキレイにエリア分け

されているのであれば、「値が1であるようなpixelの座標値をすべて平均する」と、

それも楕円の中心座標になります。

一方で、楕円が、「点列」のような形で表現されているのであれば、もうちょっと違った

アプローチが必要ですが。

（楕円の中心を求めるのに、楕円に内接する平行四辺形でなければならない理由は不明ですが、

若干面倒な気がします。やはり、その方法での解が必要でしょうか？）

[そら]

M.Shiraishi 様

すばらしいアイディアとVI を作成して頂きありがとうございました。

大変、助かりました。

今後、参考にさせて頂きます。

楕円の中心座標の求め方についてなのですが、楕円に内接する平行四辺形を使用しなくても大丈夫です。

私は、現在、瞳孔の研究に携わっているのですが、その際に、瞳孔をカメラで撮影し、撮影した瞳孔を楕円と見立てて瞳孔中心座標(楕円の中心座標)をLabVIEWを使用して求めようとしています。

今回、M.Shiraishi 様に頂いたVI を参考にしてプログラムを作成することができると思います。

ありがとうございます。

おこがましいのですが、もう一つ質問がございます。

瞳孔を撮影する際に、瞳孔にまぶたやまつげが重なり、瞳孔が半分しか撮影できない場合がございます。

私と致しましては、半分しか撮影できない場合の、瞳孔中心座標も求めたいと考えています。

求める方法と致しましては、半分の楕円から、楕円全体を推定(楕円近似)を行い、中心座標を求めようと考えています。

この方法に関しましても、プログラムを上手く作成することができていません。

目標は、半分ほど掛けた楕円において、楕円の中心座標を求めることです。

私が考えた方法以外の方法でも構いませんので、何か、アイディア、アドバイス等ご教授頂けないでしょうか？

何卒、よろしくお願いします。

（参考までに、瞳孔が半分しか撮影できなかった際の画像を添付します。）

[M.Shiraishi]

そら様

もう少し、楕円を求めるところを汎用的にしてみました。

(1) 110616-sample1.vi

(2) 110616-sample2.vi

(3) 110616-solve-ellipse-from-points-by-mse.vi　（サブvi）

---

(1) は、楕円の輪郭上にある点の座標から、最小二乗法で楕円の式 x2+Ay2+Bxy+Cx+Dy+E=0 の

パラメータ A,B,C,D,E などを求めます。

こちらは、前回、「楕円が点列で表されている場合は…」にあたりますが、楕円が部分的にしか存在しない

場合にも使えます。

そのまま実行すると、20点で楕円が描かれ、画面右側にその20点で最小二乗フィッティングしたパラメータが

表示されます。

ぴったり楕円上に乗っていますから、よほどひしゃげた楕円で無い限り、描画に使ったパラメータがそのまま

得られていると思います。

次に、「選択」スイッチを「下1/3」にして実行すると、楕円の点は下1/3だけになりますが、それでもそれらの

点だけを使ったフィッティング結果は、よほどひしゃげてない限りパラメータが再現できると思います。

（元々楕円上にあった点列であれば、かなり限られた楕円弧でも結構再現できるようです）

---

(2) は、ご要望の「欠けた瞳」画像に対して楕円をフィットしたい場合ですが、欠けた瞳の

「楕円弧」にあたる部分から点をいくつか抜き出して、点列から楕円フィットするサブviを呼び出します。

なお、欠け方によっては、輪郭のどこが元の楕円弧で、どこが不要な部分か判別するのは難しいため、

予め、「欲しい部分だけ、あるY座標でスパッと切り取った部分」を想定しています。

実際に使うときは、欠けた瞳から、不完全な部分をスパッと切り取っておいてください。

（ここでは、楕円の下部が残っている、という想定です）

viは、楕円の切り取られた部分から、Y座標を10箇所くらい選び出し、そのY断面上の輪郭の

左右の点を選び、合計20個の点で楕円フィットしています。

こちらは、元の描画は正確な楕円の式から作っていますが、輪郭の選び出しはピクセル単位なので、

フィッティング結果に多少の誤差が乗ります。

今回の例だと、最低、全体の1/4は残っていないと、さすがにかなりおかしなフィッティングになります。

ただ、ご要望の画像の瞳は半分以上残っていそうですので、それなりにご利用いただけると思います。

なお、切り取らなければ、完全な楕円に対してフィットが可能です。

---

(3) は、(1)・(2)で呼び出しているフィッティングのサブviです。

楕円の中心だけでなく、a軸半径、b軸半径、a軸のX軸からの回転角　が出てきます。

加えて、x2+Ay2+Bxy+Cx+Dy+E=0 で表した場合のA～Eも得られます。

お試し下さい。

[そら]

M.Shiraishi 様

VI を作成して頂き、また、詳しい説明、誠に有難うございます。

目的の動作達成に近づく事ができました。

今回、プログラムを拝見させて頂いたところ、疑問点が生じました。

作成して頂いたプログラムでは、Ｙ断面をＮ本選び出すのに、楕円方程式を利用しています。

これは、楕円方程式を使用しなければならないのでしょうか？

現在、あらかじめ楕円方程式を把握することができません。

そこで、楕円方程式を使用しないで、Ｙ断面をＮ本選び出す方法はありますでしょうか？

私の、理解不足かもしれませんが、よろしくお願い致します。

[M.Shiraishi]

そら様

>作成して頂いたプログラムでは、Ｙ断面をＮ本選び出すのに、楕円方程式を利用しています。

利用していません。

（利用したのでは無意味ですから・・・これから求めたいものを既知の前提にしたものでは

　役に立ちませんので）

※本命プログラムの 110616-sampl2.vi について述べています。

添付のように、Y断面をN本選んでいるところは、「部分楕円の存在するY座標範囲」を、

ある方法で直接求めています。

これは、楕円ではないところを0、楕円とその内部を1　としたとき、Y座標を順番に見ていって

（添付の赤枠内の最初のFor）、「横1列中に1があるか」（＝その横1列は楕円を含むか）の

True/Falseを取り出しますと、楕円の無いY座標はFalse、楕円のあるY座標はTrueになります。

これを、1要素ずらしてエッジ検出（ANDの片方をNOTにする）することで、「0→1になったところ」と

「1→0になったところ」が求められます。

前者は、Yをスキャンしたとき楕円が始まるY座標であり、本ソフトでは、事前に「切り取りを行ったY座標」に

相当します。後者は、楕円が終わるY座標であり、これはもともとの楕円の「最下端」を楕円方程式によらず

求めているということです。

このようにして、画像から直接「楕円の存在するY範囲」を求めましたから、その中からN本均等に

Y断面を選んでいるわけです。

Y断面を選ぶのに、（楕円を作るために使った）楕円方程式を使っているわけではないので、御確認下さい。

---

なお、サンプルの 110616-sample1.vi は、「下1/3」を作り出すのに、楕円を作った方程式を

使っていますが、これは「理論上生成した楕円から点列を作り、その点列の一部から、最小二乗法で

楕円の偏った部分だけからでも元の楕円情報を復元できますよ」という理論的な確認のために

作ったものです。

お使いになるときは、元の瞳の画像は2D配列的でしょうから、110616-sample2.vi から構築してくださいね。

[そら]

M.Shiraishi 様

質問に答えて頂き有難うございました。

説明して頂いたことについて、理解することができました。

わかりやすい説明有難うございます。

また、新たな質問があります。

作成して頂いたプログラムは、楕円でないところを0、楕円とその内部を1としています。

私は、色々と試行錯誤しているのですが、以前示した様な瞳孔画像において楕円でないところを0、楕円とその内部を1とすることができていません。

何か方法、アドバイス等ございますでしょうか？

ご教授のほど、よろしくお願い致します。

[M.Shiraishi]

今回の件に関係して、こんなものを見つけました。

http://www.vision.cs.chubu.ac.jp/04/pdf/pupil02.pdf　（添付PDF）

『高速楕円検出に基づく眼球回旋運動の計測』（中部大学、株式会社デンソー）

　第12回画像センシングシンポジウム予稿集　より

こういった感じの研究をされていますか？

このPDFのP3（予稿集のP560）に、そら様の最初の「図1」と同じ絵が、図4として描かれています。

楕円中心検出のために平行四辺形をどう使おうとしたのか、意図がだいたい分かりました。

また、最初の絵の中の平行四辺形の矢印の意味がだいたい分かりました。

本論文には、

「欠損した楕円にも対応できるように（注：ここでいう欠損とは、睫毛などで断片的に輪郭が細かく

切れている欠損ということか）、上下、左右から順に平行四辺形を求めていき、平行四辺形が求められる

場合のみ交点の画素に投票していく。最終的に投票数が最大の画素を楕円中心として算出する。」

とあり、その下に、アルゴリズム的にそのやり方が示されています。

一方で、平行四辺形法は、P6の図12に示すように、瞳の約半分が欠けたような場合は使えない、とありますね。

最小二乗フィットならうまく輪郭点が5点以上選べれば（だいたい10点以上あれば）、半分欠けてても

1/3くらい残っているだけでもフィッティングできるのですが。

P4の3.1.4で、式(3)は楕円のパラメータとしては1つ過剰ですね。

（全体をAで割ればよく、x^2の係数は1でよい）

・・・それで、平行四辺形法にしろ最小二乗法にしろ、瞳領域を割り出すために、ご質問の

「瞳の領域の検出」が必要になりますが、上記論文のP3にも

「楕円推定を行うには画像中から楕円の輪郭候補点を抽出しなければならない。

近赤外画像では、瞳孔領域が黒く浮き出るため二値化処理によって瞳孔領域を検出する。」

とサラっと書いている部分を、具体的に、どう実装したら良いのか？ということですよね。

二値化も、うまくいくようで、現実には画像全体に明るさのグラデーションがかぶって、単一の閾値では

瞳の片側はいいけど反対側は白飛びするとか、逆に片側エッジを捉えると反対側は消えてなくなるとか、

なかなかうまくいかないことが多そうです。

瞳は画像の中で最も暗いエリアなので、うまくコントラスト微分値なんかを使って検出したいところですが、

むしろ画像処理的な専門家にアドバイスを仰ぎたい領域です。

[そら]

M.Shiraishi 様

私は、提示して頂いた論文を含め、その他多くの論文を参考に、独自の視線検出システムを開発しようと研究を行っております。

楕円に関する質問はその為のものです。

今回は、たくさんの質問に親切かつ分かりやすく答えて頂きありがとうございました。

さらに、VI作成、瞳孔についての論文まで調査して頂き、心から感謝しております。

また、最後の質問では、Lab VIEW関連以外の質問をしてしまい申し訳ございません。

今後、M.Shiraishi 様から頂いたアドバイス、アイディアを元に、目標の動作達成に向けて励んでいきたいと思います。

まだまだ、Lab VIEWに関して未熟なので、機会があれば今後ともよろしくお願い致します。