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

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

行列Mを

として、ベクトル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;
}

解は

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

じゃないでしょうか？？

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

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

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

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

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