3次元 円の中心座標の求め方

Question

前提・実現したいこと

C言語で3次元の3点(x1,y1,z1),(x2,y2,z2),(x3,y3,z3)を通る
円の中心座標を求めたいです。

発生している問題・エラーメッセージ

歯食べたところ求められそうな連立方程式の式は出てきたのですが、
それをC言語で実現する方法がわかりません。
また内積や外積等を使うとかんたんにできるとも見ましたが、
よくわかっていません。

簡単なサンプルがあればそれを動かしながら理解していけると思っているのですが、見当たらずこちらで質問させていただきました。

利用環境はgccなどC言語であれば何でもOKです。

よろしくお願いします。

Answer 1

中心P(x, y, z) の座標を求めたいのなら、
連立一次方程式を作って解けばいいだけの話でしょう。

#include <stdio.h>  // printf
#include <math.h>   // fabs

void gauss(double a[3][3 + 1])  // ガウスの消去法
{
	for (int i = 0; i < 3; i++) {
		double m = 0, t;
		int pivot = i;
		for (int j = i; j < 3; j++)
			if (fabs(a[j][i]) > m)
				m = fabs(a[j][i]), pivot = j;
		if (pivot != i)
			for (int j = 0; j <= 3; j++)
				t = a[i][j], a[i][j] = a[pivot][j], a[pivot][j] = t;
	}
	for (int k = 0; k < 3; k++) {
		double p = a[k][k];
		a[k][k] = 1;
		for (int j = k + 1; j <= 3; j++)
			a[k][j] /= p;
		for (int i = k + 1; i < 3; i++) {
			double q = a[i][k];
			for (int j = k + 1; j <= 3; j++)
				a[i][j] -= q * a[k][j];
			a[i][k] = 0;
		}
	}
	for (int i = 3; --i >= 0; )
		for (int j = 3; --j > i; )
			a[i][3] -= a[i][j] * a[j][3];
}

typedef struct Point {
	double x, y, z;
} Point;

Point getCenter(Point A, Point B, Point C)
{
	double M[3][3 + 1];  // Matrix (方程式の係数)

	// P(x, y, z) は中心の座標
	// |AP| = |BP|
	// (x - A.x)^2 + (y - A.y)^2 + (z - A.z)^2
	//    = (x - B.x)^2 + (y - B.y)^2 + (z - B.z)^2
	M[0][0] = 2*(B.x - A.x);
	M[0][1] = 2*(B.y - A.y);
	M[0][2] = 2*(B.z - A.z);
	M[0][3] = B.x*B.x + B.y*B.y + B.z*B.z - A.x*A.x - A.y*A.y - A.z*A.z;

	// |AP| = |CP|
	// (x - A.x)^2 + (y - A.y)^2 + (z - A.z)^2
	//    = (x - C.x)^2 + (y - C.y)^2 + (z - C.z)^2
	M[1][0] = 2*(C.x - A.x);
	M[1][1] = 2*(C.y - A.y);
	M[1][2] = 2*(C.z - A.z);
	M[1][3] = C.x*C.x + C.y*C.y + C.z*C.z - A.x*A.x - A.y*A.y - A.z*A.z;

	double AB[3] = { B.x - A.x, B.y - A.y, B.z - A.z };  // ベクトルAB
	double AC[3] = { C.x - A.x, C.y - A.y, C.z - A.z };  // ベクトルAC
	double ABxAC[3] = {   // 外積ABxAC は平面ABCの法線ベクトル
		AB[1]*AC[2] - AB[2]*AC[1],
		AB[2]*AC[0] - AB[0]*AC[2],
		AB[0]*AC[1] - AB[1]*AC[0],
	};
	// P は 平面ABC上、法線ベクトルABxAC とベクトルAP は直交、内積が 0
	// ABxAC[0](x - A.x) + ABxAC[1](y - A.y) + ABxAC[2](z - A.z) = 0
	M[2][0] = ABxAC[0];
	M[2][1] = ABxAC[1];
	M[2][2] = ABxAC[2];
	M[2][3] = ABxAC[0]*A.x + ABxAC[1]*A.y + ABxAC[2]*A.z;

	gauss(M);

	Point P = { M[0][3], M[1][3], M[2][3] };
	return P;	
}

int main(void)
{
	Point A = { 2, 2, 5 }; // Point A = { x1, y1, z1 };
	Point B = { 1, 1, 5 }; // Point B = { x2, y2, z2 };
	Point C = { 3, 1, 5 }; // Point C = { x3, y3, z3 };

	Point P = getCenter(A, B, C);
	printf("(%g, %g, %g)\n", P.x, P.y, P.z);
}

Comments

[santyago]

サンプルありがとうございました。
連立方程式をC言語で解く方法がわかってませんでした。
こちらを参考に理解していきたいと思います。

[kazuma-s]

> 歯食べたところ求められそうな連立方程式の式は出てきたのですが、

その調べた連立方程式というのを、質問に追記してもらえませんか？
質問は編集できます。

Answer 2

（３点が一直線に並んでいなければ，という条件が前提としてあるでしょうが）
３点を頂点とした三角形の外心（外接円の中心）を求めればよいのではないでしょうか．

外接円の求め方はググれば見つかるかと思います．

Comments

[fana]

質問文の雰囲気から，ガチで連立方程式を扱うよりもマシな方法を求めているのかと思ったのだけど，ちがったみたい．

Answer 3

歯食べたところ求められそうな連立方程式の式は出てきたのですが、

それをC言語で実現する方法がわかりません。

まずは連立方程式を数学的に解いて、一般解を求めましょう。そうすればあとはその解をC言語で書くだけです。