多点の近くを通る円の中心


ある任意の4点を通る円の曲率を求めようとしています.
4点をきちんと通る円であることはまれでしょうから,最小二乗法を用いて
4点のできるだけ近くを通る円を考えています.問題設定はここと同じです.
http://oshiete.eibi.co.jp/kotaeru.php3?q=1110779

それで,上記URLの回答No.3 No.6の回答を元に計算していたのですが,
No.6の中盤で

> X^2*A + XY*B + X*C + X(X^2+Y^2) = 0
> XY*A + Y^2*B + Y*C + Y(X^2+Y^2) = 0
> X*A + Y*B + C + (X^2+Y^2) = 0
> という式が得られます。これをA,B,Cを変数とする連立方程式として解いて

という部分を実行しようとしてもここは,どうやっても解けない気がします.
それで,この部分の解決法,あるいはまったく違う視点から,4点のできるだけ
近くを通る円を求めるやり方をご存知でしたら,教えていただきたく思います.

よろしくお願いいたします.

回答の条件
  • 1人3回まで
  • 登録:2008/01/21 03:57:51
  • 終了:2008/01/28 04:00:03

回答(3件)

id:ita No.1

ita回答回数203ベストアンサー獲得回数472008/01/21 08:43:41

ポイント27pt

未知数A,B,Cで指揮が3つあるから解けると思いますが。

行列Mを

X^2 XY X
XY Y^2 Y
X Y 1

として、ベクトルvを(X(X^2+Y^2), Y(X^2+Y^2), (X^2+Y^2))とおき、

M (A,B,C)=-v を解きます。具体的にはMの逆行列M^-1を-v に掛けます。3x3の逆行列を計算するコードを以下に示します。

#include <math.h>

typedef double Vec3[3];
typedef double Mat3[3][3];

extern double mat_det3(Mat3 m)
{
  double m11=m[0][0];
  double m22=m[1][1];
  double m33=m[2][2];

  double m12=m[0][1];
  double m21=m[1][0];
  double m13=m[0][2];
  double m31=m[2][0];
  double m23=m[1][2];
  double m32=m[2][1];
  return m11*(m22*m33-m23*m32)-m12*(m21*m33-m23*m31)+m13*(m21*m32-m22*m31);
}

extern int mat_inv3(Mat3 m, Mat3 inv)
{
  double m11,m22,m33,m12,m13,m23, m21, m31,m32;
  double det;
  int i,j;

  m11=m[0][0];
  m22=m[1][1];
  m33=m[2][2];
  m12=m[0][1];
  m21=m[1][0];
  m13=m[0][2];
  m31=m[2][0];
  m23=m[1][2];
  m32=m[2][1];

  det= mat_det3(m);

  inv[0][0]= m22*m33 -m23*m32;
  inv[1][1]= m11*m33 -m13*m31;
  inv[2][2]= m11*m22 -m12*m21;

  inv[0][1]= m13*m32 - m12*m33;
  inv[1][0]= m31*m23 - m21*m33;
  inv[0][2]= m12*m23 - m13*m22;
  inv[2][0]= m21*m32 - m31*m22;
  inv[1][2]= m13*m21 - m11*m23;
  inv[2][1]= m31*m12 - m11*m32;

  det = 1.0 / det;

  for(i=0;i<3;i++)
  for(j=0;j<3;j++) inv[i][j] *=det;
  return 0;
}
id:violin

ita様


ご回答ありがとうございます.


小生,大学以上の数学はもとよりプログラミングも入門者以前なので,フリーの数式処理ソフトMaximaやMupadを使って実行いたしました.ita様のご回答のように逆行列を用いて実行しましたがエラーとなりました.rankを求めると1しかないことが原因でしょうか.


問題の式を変形させるとこうなりますよね.

X(X*A+Y*B+C+(X^2+Y^2))=0

Y(X*A+Y*B+C+(X^2+Y^2))=0

(X*A+Y*B+C+(X^2+Y^2))=0


ちなみに,式を普通に求めようとするとこうなります.

http://picasaweb.google.com/mit.sfc/HaSFzI/photo#515800706731316...

2008/01/22 17:56:15
id:Rec No.2

Rec回答回数10ベストアンサー獲得回数32008/01/21 17:28:40

ポイント27pt

解は

A=0, B=0, C=-X^2-Y^2

じゃないでしょうか??

上の連立方程式に代入してみると全て満たします。

id:violin

Rec様


ご回答ありがとうございます.


A,Bを0としたとき,3つの式に代入すると

C=-(X^2+Y^2)/X

C=-(X^2+Y^2)/Y

C=-(X^2+Y^2)

となり,Cは不定だと思いますが.

これは私の知識不足で,定まるというのでしょうか.


今取り組んでいる方式に限らず,4点を近似して通る円の曲率を求める式がありましたら,そちらを教えていただきたく思います.どうぞよろしくお願いいたします.

2008/01/22 17:56:52
id:ita No.3

ita回答回数203ベストアンサー獲得回数472008/01/22 14:58:20

ポイント26pt

しまった、行列のランクは明らかに1ですね( (1,X,Y)×(1,X,Y)と書ける)。

SF作家のグレッグ・イーガンがこの問題に関する記事を書いてました。

http://groups.google.com/group/comp.graphics.algorithms/browse_t...

x, y座標の分散、共分散を求め、分散・共分散行列(2x2)の逆行列を計算するとかで計算できるようです。

id:violin

ita様


やはりrankが1だからですか。

興味深いURIもありがとうございます。時間ができたときに拝見させていただきます。


この問題はこの方法では解けそうもないですね。

ここでの質問者や回答者は果たして解けたのか、ちょっと気になります。。。

http://oshiete.eibi.co.jp/kotaeru.php3?q=1110779


この問題に私自身が時間を割ける余裕がだんだんなくなってきましたので、質問文のような円を求める式をご存じでしたら、絶対にこれで出ますよということをばからしいくらい丁寧にお答えくだされば幸いです。


みなさま、よろしくお願い申し上げます。

2008/01/22 18:03:13
  • id:splint
    http://picasaweb.google.com/mit.sfc/HaSFzI/photo#5158007067313166674

    これが出ているならABCも出るのでは? XもYもXYも求まるわけですし。
  • id:ita
    あ、やっぱり1の方法で解けます。
    X^2 とかはx二乗の和という意味なんで因数分解してはだめです。
    数式処理するときはsxxとか書くといいかも。
    solve( sxx*a+sxy*b+x*c=sx*(sxx+syy),
    sxy*a+syy*b+sy*c=sy*(sxx+syy),
    sx*a+sy*b+c=sxx+syy, (a,b,c))

    でどうでしょう。

この質問への反応(ブックマークコメント)

「あの人に答えてほしい」「この質問はあの人が答えられそう」というときに、回答リクエストを送ってみてましょう。

これ以上回答リクエストを送信することはできません。制限について

絞り込み :
はてなココの「ともだち」を表示します。
回答リクエストを送信したユーザーはいません