Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this community
Can someone please help in finding an analytic expression for the center of a circle through 3 given points by modifying ( Solve in place of NSolve ) a conventional numerical approach above? The program hangs and some more coding is needed. Or, from the following we may need to work with just 3 input points in 3-space.
Thanks in advance. Shall provide your answer link in my referenced answer above.
We first solve the problem for 2D and subsequently for 3D.
For 2D we define the following function that calculates the center of 3 points. For this we first calculate the mid points between the original points. Then we define the perpendicular bisectors of the lines as functions of some "li". The crossing point of the bisectors is the searched for center.
Clear[getCenter]
getCenter[pts_] := Module[{ms, pb, sol, cen, l1, l2, l3},
ms = (pts + RotateLeft[pts])/2; (*mid points of the sides*)
pb = ms + {l1, l2, l3} ({{0, -1}, {1, 0}} . # & /@ (pts - RotateLeft[pts])); (* perpendicular bisectors*)
sol = Solve[pb[[1]] == pb[[2]], {l1, l2}]; (*li of crossing point*)
cen = pb[[1]] /. sol[[1]] (*center*)
]
Now for 3D. The idea is to rotate the 3 points in a plane parallel to the x/y plane. If we then eliminate the z-coordinate we have a 2D problem. Toward this aim, we first rotate calculate the normal to the plane through the 3 points: "no". The searched for rotation is then defined by rotating the normal so that is parallel to the z-axis. The rotated points will all lie in a plane parallel to x/y plane. To get the center of the circle in this plane we use the function from above. To rotate the center back to the original points we need the inverse rotation, what is the same as the transposed as the rotation is a orthonormal transformation.
Here is an example:
ps = {p1, p2, p3} = RandomReal[{-1, 1}, {3, 3}]; (*3 points*)
no = Cross[p2 - p1, p3 - p1]; (*normal to plane*)
rot = RotationTransform[{no, {0, 0, 1}}][[1, ;; 3, ;; 3]]; (*rotation matrix*)
rps = rot . # & /@ ps; (*rotated points*)
cen = Transpose[rot] . Append[getCenter[rps[[All, ;; 2]]], rps[[1, 3]]] (*center of circle through points*)
Graphics3D[{PointSize[0.02], Point[ps], Green, Point[cen], Blue,
Opacity[0.1], Cylinder[{cen, cen + 0.001 no}, Norm[ps[[1]] - cen]]},
Axes -> True, AxesLabel -> {"x", "y", "z"}]
You can get an expression for the centre in terms of 3D positions of 3 points on the target. The centre is defined in terms of 3 planes:
Put the coordinates in a 3x3 matrix
M1 = Through /@ Array[{X, Y, Z}, {3}]
(* {{X[1], Y[1], Z[1]}, {X[2], Y[2], Z[2]}, {X[3], Y[3], Z[3]}} *)
We can define a test for (x,y,z) co-planar in terms of projective coordinates
Append[M1, {x, y, z}]
(* {{X[1], Y[1], Z[1]}, {X[2], Y[2], Z[2]}, {X[3], Y[3], Z[3]}, {x, y,
z}} *)
M = Append[#, 1] & /@ %
(* {{X[1], Y[1], Z[1], 1}, {X[2], Y[2], Z[2], 1}, {X[3], Y[3], Z[3],
1}, {x, y, z, 1}} *)
planar = FullSimplify[Det[M] == 0]
(* z (X[3] (Y[1] - Y[2]) + X[1] (Y[2] - Y[3]) + X[2] (-Y[1] + Y[3])) +
X[3] Y[2] Z[1] + x Y[3] Z[1] + x Y[1] Z[2] +
X[1] Y[3] Z[2] + (X[2] Y[1] - X[1] Y[2] + x (-Y[1] + Y[2])) Z[3] +
y (X[3] (-Z[1] + Z[2]) + X[2] (Z[1] - Z[3]) +
X[1] (-Z[2] + Z[3])) ==
x Y[2] Z[1] + X[2] Y[3] Z[1] + X[3] Y[1] Z[2] + x Y[3] Z[2] *)
Now find the condition that (x,y,z) is equidistant from the 3 points
# - {x, y, z} & /@ M1
(* {{-x + X[1], -y + Y[1], -z + Z[1]}, {-x + X[2], -y + Y[2], -z +
Z[2]}, {-x + X[3], -y + Y[3], -z + Z[3]}} *)
Equal @@ (# . # & /@ %)
(* (-x + X[1])^2 + (-y + Y[1])^2 + (-z + Z[1])^2 == (-x + X[2])^2 + (-y +
Y[2])^2 + (-z + Z[2])^2 == (-x + X[3])^2 + (-y + Y[3])^2 + (-z +
Z[3])^2 *)
equidistant = Simplify[Expand /@ %]
(* 2 x X[1] + X[2]^2 + 2 y Y[1] + Y[2]^2 + 2 z Z[1] + Z[2]^2 ==
X[1]^2 + 2 x X[2] + Y[1]^2 + 2 y Y[2] + Z[1]^2 + 2 z Z[2] &&
2 x X[2] + X[3]^2 + 2 y Y[2] + Y[3]^2 + 2 z Z[2] + Z[3]^2 ==
X[2]^2 + 2 x X[3] + Y[2]^2 + 2 y Y[3] + Z[2]^2 + 2 z Z[3] *)
Mathematica has no difficulty in solving this (but the result is quite large). Note that much of the complexity comes because symbolic solution of 3x3 matrices gives complicated expressions (in general)
solution1 = Solve[equidistant && planar, {x, y, z}];
To verify the solution, define a circle
centre = {3, 4, 5};
radius = 2;
normal = {3/13, 4/13, 12/13};
normal . normal == 1
(* True *)
Get Mathematica to generate 3 points on the circle
pts = Block[{v = {x, y, z}},
FindInstance[
normal . v == normal . centre && (v - centre) . (v - centre) ==
radius^2, v, Reals, 3]
];
Thread /@ Thread[M1 -> ({x, y, z} /. pts)] // Flatten
(* {X[1] -> 33/26, Y[1] -> (4295 - 39 Sqrt[535])/1040,
Z[1] -> (5605 + 13 Sqrt[535])/1040, X[2] -> 119/26,
Y[2] -> (4037 - 39 Sqrt[879])/1040,
Z[2] -> (4831 + 13 Sqrt[879])/1040, X[3] -> 5/2,
Y[3] -> 1/80 (323 - 3 Sqrt[2391]), Z[3] -> 1/80 (409 + Sqrt[2391])} *)
Verify that this locates the solution correctly.
(solution1 /. %) // Simplify
(* {{x -> 3, y -> 4, z -> 5}} *)
Method 1
eqn = {
(x - x1)^2 + (y - y1)^2 + (z - z1)^2 ==
(x - x2)^2 + (y - y2)^2 + (z - z2)^2 ==
(x - x3)^2 + (y - y3)^2 + (z - z3)^2 ,
Det[{
{x, y, z, 1},
{x1, y1, z1, 1},
{x2, y2, z2, 1},
{x3, y3, z3, 1}
}] == 0
};
sol = Solve[eqn, {x, y, z}] // Simplify; // AbsoluteTiming
LeafCount[center1 = {x, y, z} /. sol[[1]]]
{1.27057, Null}
3315
Method 2
{a1,b1,c1,d1}={-2*x1+2*x2,-2*y1+2*y2,-2*z1+2*z2,x1^2-x2^2+y1^2-y2^2+z1^2-z2^2};
{a2,b2,c2,d2}={-2*x1+2*x3,-2*y1+2*y3,-2*z1+2*z3,x1^2-x3^2+y1^2-y3^2+z1^2-z3^2};
{a3,b3,c3,d3}=Det/@{{{1,1,1},{y1,y2,y3},{z1,z2,z3}},{{x1,x2,x3},{1,1,1},{z1,z2,z3}},{{x1,x2,x3},{y1,y2,y3},{1,1,1}},-{{x1,x2,x3},{y1,y2,y3},{z1,z2,z3}}};
{d,dx,dy,dz}=Det/@{{{a1,b1,c1},{a2,b2,c2},{a3,b3,c3}},{{-d1,b1,c1},{-d2,b2,c2},{-d3,b3,c3}},{{a1,-d1,c1},{a2,-d2,c2},{a3,-d3,c3}},{{a1,b1,-d1},{a2,b2,-d2},{a3,b3,-d3}}};
LeafCount[center2={dx,dy,dz}/d]
center1 - center2 // Simplify
927
{0,0,0}
