1
$\begingroup$

I'd like to find the best fit coefficients "a" and "c" for a toroidal surface, to a list of {x,y,z} points which are approximately on the surface of the torus. I believe I need LinearModelFit for this, the function for the torus: z^2 == a^2 - (c - (x^2 + y^2)^(1/2))^2 and the data:

data={{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},{x5,y5,z5}};   (* these would be passed as actual xyz values not variables. My fit function is more complicated than below so the actual values I have are not appropriate *)

In which case the command should look something like this but I'm clearly not expressing the torus function as a model correctly for LinearModelFit. Also I will need to include some constraints, but once I have the syntax and functional form for the model correct I should be able to manage the constraints and starting values:

LinearModelFit[data, z^2 == a^2 - (c - (x^2 + y^2)^(1/2))^2,{a,c,,x,y,z}]

This answer How to fit a surface to 3D data in Mathematica? is helpful but I couldn't see how to express the model for a torus in an equivalent form.

$\endgroup$
2
  • 1
    $\begingroup$ If your torus is rotated or translated to a particular offset, then you could make use of a KLDecomposition of the points to set them at the right orientation and position them in the center. See this answer mathematica.stackexchange.com/a/225226/72682. If this works with minimal noise you would only need to focus on fitting the two radii. The resulting torus could then be rotated and translated back using the KL matrix. Also, there's a torus here in case you want to visualize it: Torus = ResourceFunction["Torus"]; resources.wolframcloud.com/FunctionRepository/resources/Torus $\endgroup$
    – flinty
    Aug 5 '20 at 1:21
  • $\begingroup$ Torus fitting appears here on page 638 if you're interested. $\endgroup$
    – flinty
    Aug 5 '20 at 2:09
3
$\begingroup$

Assuming that data is obtained from translated and/or rotated models

Clear[a, c]
p = {x, y, z};
P = {X, Y, Z};
Txyz = z^2 - a^2 + (c - (x^2 + y^2)^(1/2))^2 /. Thread[p -> RollPitchYawMatrix[{alpha, beta, gamma}].(P - {x0, y0, z0})];

now using the excellent script from @flinty for random data generation

SeedRandom[1];
Torus = ResourceFunction["Torus"];
testTorus = Torus[{4, -2, 6}, {19, 4}];
(*pts on a torus plus some noise*)

pts = RandomPoint[DiscretizeGraphics@testTorus, 300] + RandomVariate[NormalDistribution[0, .5], {300, 3}];
gr1 = Graphics3D[Point@pts, Axes -> True];

and then following with a minimization procedure

error = Sum[(Txyz /. Thread[P -> pts[[k]]])^2, {k, 1, Length[pts]}];
sol = NMinimize[{error, -Pi <= alpha <= Pi, -Pi <= beta <= Pi, -Pi <= gamma <= Pi}, {a, c, x0, y0, z0, alpha, beta, gamma}]

and the results

gr0 = ContourPlot3D[(Txyz /. sol[[2]]) == 0, {X, -20, 20}, {Y, -20, 20}, {Z, -20, 20}, ContourStyle -> {Yellow, Opacity[0.2]}, Mesh -> None, BoundaryStyle -> None];
Show[gr1, gr0]

enter image description here

$\endgroup$
3
  • $\begingroup$ The use of NMinimize for this is much appreciated as it also highlights how the fit is obtained. Also adaptable to many other problems. $\endgroup$
    – DrBubbles
    Aug 5 '20 at 17:48
  • $\begingroup$ Very nice. But I don't see where "b" enters the "error" and so shouldn't it be omitted from the list of variables in NMinimize? $\endgroup$
    – mef
    Aug 6 '20 at 19:21
  • $\begingroup$ @mef Yes. The variable b should not participate in the modeling. It is now corrected. Thanks. $\endgroup$
    – Cesareo
    Aug 6 '20 at 20:13
2
$\begingroup$

Shape fitting is quite difficult and has robustness issues without taking special care, but you can get quite far without any fitting. If your torus data is not rotated and oriented directly up, then the BoundingRegion commands provide very good parameters:

SeedRandom[1];
Torus = ResourceFunction["Torus"];
testTorus = Torus[{4, -2, 6}, {19, 4}];
(* pts on a torus plus some noise *)
pts = RandomPoint[DiscretizeGraphics@testTorus, 300] + 
   RandomVariate[NormalDistribution[0, .5], {300, 3}];

minball = BoundingRegion[pts, "MinBall"];
minbox = BoundingRegion[pts, "MinCuboid"];
pos = minball[[1]];
radius = minball[[2]];
holeradius = Min[EuclideanDistance[pos, #] & /@ pts];
piperadius = Min[Abs[minbox[[1]] - minbox[[2]]]]/2;
Graphics3D[{Point@pts, Opacity[.5], 
  Torus[pos, {(radius + holeradius), piperadius}]}]

torus fit

If you torus data is oriented, then you should use the KarhunenLoeveDecomposition to get it into a face-up and centered form first, then fit, then rotate and translate back. Let me know if this is desired too.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.