# Is there a method faster than FindFit to fit conical elliptical cylinder to point data?

I have a code that fits conical elliptic cylinders to point data. It utilizes FindFit.

Usually, I need to fit over 1,500 "cylinders"**, each are represented by ~ 100 points. The full algo (which is not presented here) makes up to 10 iterations,eliminating some "bad" points along the way, which comes up to 15,000 calls of FindFit, each takes about 0.05 seconds. This leads to over 10 minutes of calculation.

Is there a faster way, other than FindFit? Like, an order of magnitude faster..?

** These are not real cylinders. They have elliptic cross section and they tend to be thicker at one end, which is why I call them "conical elliptic cylinders".

Have I tried alternatives, myself? Yes. There is RegionFit. It runs pretty fast. But it only fitst "Cylinders" and "Cones". It does not have "Conical Elliptic Cylinder" option.

Toy model CSV Data: one "cylinder".

Point Data

Current code:

SetDirectory[NotebookDirectory[]];
points = Import["points.csv"];

(*Constraints*)
RsqMin = 0.2;
RsqMax = 0.3;
rotationMinMax = Pi/18;
translationMinMax = 0.5;
zConeCoefMinMax = 0.1;

fitConicalEllipticCylinder[pts_, RsqMin_, RsqMax_, rotationMinMax_, translationMinMax_, zConeCoefMinMax_] := Module[{fitPts, model, coeffs, rX, rY, rZ, t, transformation, rZ1, rXY, tZ, x, y, z},
rX = RotationTransform[angX, {1, 0, 0}];
rY = RotationTransform[angY, {0, 1, 0}];
rZ = RotationTransform[angZ, {0, 0, 1}];
t = TranslationTransform[{tX, tY, 0}];
transformation = Composition[rZ, rY, rX, t];
model = (transformation[{x, y, z}][])^2/Asq + (transformation[{x, y, z}][])^2/Bsq - (1 - transformation[{x, y, z}][]*zC)^2;
coeffs = {Asq, Bsq, angX, angY, angZ, tX, tY, zC};
fitPts = {#1, #2, #3, 0} & @@@ pts;
FindFit[fitPts, {model, {RsqMin < Asq < RsqMax, RsqMin < Bsq < RsqMax, -rotationMinMax < angX < rotationMinMax, -rotationMinMax < angY < rotationMinMax, -rotationMinMax < angZ < rotationMinMax, -translationMinMax < tX < translationMinMax, -translationMinMax < tY < translationMinMax, -zConeCoefMinMax < zC < zConeCoefMinMax}}, coeffs, {x, y, z}]
]

fitResult = fitConicalEllipticCylinder[points, RsqMin, RsqMax, rotationMinMax, translationMinMax, zConeCoefMinMax] This took ~ 0.06 seconds on my machine.

a = Sqrt[fitResult[[1, 2]]];
b = Sqrt[fitResult[[2, 2]]];
coneC = fitResult[[8, 2]];

pointplog = ListPointPlot3D[points, Axes -> True, BoxRatios -> Automatic];
parametricplot = ParametricPlot3D[{Cos[u]*a*(1 - v*coneC), Sin[u]*b*(1 - v*coneC), v}, {u, 0, 2 Pi}, {v, Min[points[[All, 3]]], Max[points[[All, 3]]]}, Mesh -> None, PlotStyle -> {FaceForm[LightBlue], Opacity[0.5]}];

Show[pointplog, parametricplot] • Right conical elliptic cylinders or oblique conical elliptic cylinders? Aug 19 at 11:30
• It’s scan data, so I don’t really know. I treat them as RIGHT. Aug 19 at 11:55
• There is BoundingRegion but it doesn't seem to have an elliptic cylinder as a target for fitting. Perhaps there is a workaround or an undocumented option.
– Syed
Aug 19 at 12:05
• True. That's why I had to go with FindFit. Now, I hope to maybe make it faster. Aug 19 at 13:43
• Because angX, angY, angZ, tx, tY, and zC seem to be such small numbers (at least in this particular dataset), putting in better starting values than the default value of 1 reduces the time from around 0.06 to 0.05: coeffs = {{Asq, 0.25}, {Bsq, 0.25}, {angX, 0}, {angY, 0}, {angZ, 0}, {tX, 0}, {tY, 0}, {zC, 0}};. Not lots of savings but these starting values are within the constraints which the default values of 1 are not.
– JimB
Aug 19 at 15:31

A general quadric is the zero-set of the quadratic polynomial in $$x \in \mathbb{R}^3$$ of the form $$\tfrac{1}{2} x^T \varSigma \, x + \beta^T x = 1,$$ where $$\varSigma$$ is some symmetric $$3 \times 3$$ matrix and $$\beta \in \mathbb{R}^3$$ is some vector.

An elliptic cone is a special instance of this for which exactly one of the the three eigenvalues of $$\varSigma$$ is equal to zero smaller or equal to zero and the other two are positive. If this eigenvalue is zero then you actually have an elliptic cylinder (or some degenerate object).

So if the data points are "sufficiently clean", then it should be possible to fit a quadric to the datasets and two of the eigenvalues of $$\varSigma$$ will be positive and one will be negative or close to zero. Since this is a fitting problem, I am going to use FindFit, too. But I hope that I can just ignore to set any further constraints. That should make it a lot easier for FindFit to find a solution, as it is just an unconstrained polynomial fit. Importing your dataset into the variable pts, my code looks like this:

\[CapitalSigma] = Array[\[Sigma], {3, 3}];
\[CapitalSigma] = UpperTriangularize[\[CapitalSigma]] + Transpose[UpperTriangularize[\[CapitalSigma], 1]];
\[Beta] = Array[b, {3}];
X = Array[x, {3}];
vars = Join[DeleteCases[Flatten[UpperTriangularize[\[CapitalSigma]]],0], \[Beta]];

f = x |-> 1/2 x.\[CapitalSigma].x + \[Beta].x;
result = FindFit[Join[pts, ConstantArray[1., {Length[pts], 1}], 2], f[X], vars, X]; // RepeatedTiming // First


0.0000905272

That is quite a bit faster, isn't it? But of course, we have not achieved everything OP asked for. And we still have to see how good the solution is. So let's plot the level set:

Show[
Graphics3D[Sphere[pts, 0.02]],
ContourPlot3D[
Evaluate[f[X] == 1 /. result], {x, -1.5, 1.5}, {x, -1.5, 1.5}, {x, -1.5, 1.5}, PlotPoints -> 60]
] This looks pretty good. But the data set was also pretty clean. I don't know how well this will work with real life data.

If you look very closely on this you will see that the quadric is actually an ellipsoid with just a very long axis. This reflects the fact that one eigenvalue is positive, but quite small compared to the others:

Eigenvalues[\[CapitalSigma] /. result]


{7.98575, 7.65485, 0.0406362}

We can repair this by just killing this smallest eigenvalue (smallest by absolute value) clipping the smallest eigenvalue of $$\varSigma$$ to the interval $$]-\infty,0]$$.

{\[Lambda], U} = Eigensystem[\[CapitalSigma] /. result];
\[CapitalSigma]new = U\[Transpose] . ({\[Lambda][], \[Lambda][], Clip[\[Lambda][], {-\[Infinity], 0.}]} U);

Show[
Graphics3D[Sphere[pts, 0.02]],
ContourPlot3D[
Evaluate[
1/2 X . \[CapitalSigma]new . X + \[Beta] . X == 1 /. result], {x, -1.5, 1.5}, {x, -1.5, 1.5}, {x, -1.5, 1.5},
PlotPoints -> 60]
] We don't solve an optimization problem here, so there might be better instances of $$\varSigma$$ with exactly two nonzero eigenvalues. But this strategy is inexpensive and might perform sufficiently well.

The axis of the cone goes along U[]. So if you want to find planes to truncate this to a elliptical conical cylinder, I would look for the two tightest planes orthogonal to U[].

{min, max} = MinMax[pts.U[]];

Show[
Graphics3D[Sphere[pts, 0.02]],
ContourPlot3D[
Evaluate[
1/2 {x, y, z}.\[CapitalSigma]new.{x, y, z} + \[Beta].{x, y, z} == 1 /. result],
{x, -1.5, 1.5}, {y, -1.5, 1.5}, {z, -1.5,1.5},
PlotPoints -> 60,
RegionFunction -> Function[{x, y, z}, Evaluate[min <= {x, y, z} . U[] && {x, y, z} . U[] <= max]],
RegionBoundaryStyle -> None
]
] Comments and further directions:

Regarding OP's question on the geometric parameters from the comments: I thought a bit about it and then observed that some of my original statements were wrong. (Which does not mean that the method does not work.) The normal form of a cone is

$$x^2/a^2 +y^2/b^2 - z^2/c^2 = 0.$$

So $$\varSigma$$ can have two positive eigenvalues and a negative one, too. But the generic shape of such signature is a hyperboloid with normal form

$$x^2/a^2 +y^2/b^2 - z^2/c^2 = d$$

with $$d \neq 0$$. So I guess, the best strategy is to bring the fitted quadric in normal form (the eigenbasis of $$\varSigma$$ tells you how to rotate it and the eigenvalues are $$1/a^2$$, $$1/b^2$$, and $$-1/c^2$$. What remains to do is to complete the square to absorb the (rotated) $$\beta$$. That is you can rewrite your quadric as

$$(x-\mu)^T \, U^T \varSigma \, U \, (x- \mu) = d$$

with the orthogonal matrix $$U$$ built from the eigenvectors. Then $$U^T \varSigma \, U$$ becomes diagonal:

$$U^T \varSigma \, U = \begin{pmatrix}\lambda_1&0&0\\0&\lambda_2&0\\0&0&\lambda_3\end{pmatrix}.$$

Let's suppose that $$\lambda_1 \geq \lambda_2 > \lambda_3$$, then the eigenvalue to $$\lambda_3$$ should be the axis of the cone or cylinder. If $$\lambda_3 = 0$$, then you get an elliptic cylinder. If $$\lambda_3 < 0$$, then you get a hyperbolid if $$d \neq 0$$ and an elliptic cylinder for $$d = 0$$.

Hence I suggest to check whether $$\lambda_3$$ is positive; if yes, then set it to $$0$$. If $$\lambda_3$$ is negative then clip $$d$$ to $$0$$ to get the conical case.

The other geometric data can be extracted from the normal form. In the conical case, $$\mu$$ should be the apex of the cone. Radii/half axes and what you call "center" can be off from the normal form once you have computed {min, max} = MinMax[pts.U[]];. Basically you have to substitute the values min and max (or some translation of them due to the completion of the square) as $$z$$-values into the normal form. Then you get the equations for the cross sectional ellipses.

• Henrik, maybe there's a typo somewhere or a missing variable? I plugged my data into your code and the output says "FindFit::fitc: Number of coordinates (3) is not equal to the number of variables (1)." three times and end with "General::stop: Further output of FindFit::fitc will be suppressed during this calculation." Aug 20 at 7:56
• Oh I see. Being afk, but I think X = Array[x,3] should help. Have to look into this tonight. Aug 20 at 8:47
• Okay, I checked it in the meantime. Seems to work know. Aug 20 at 11:13
• Yes, it does! I'm impressed. It's over 200 times faster then my "OP" code, yet it catches all the rotation & translation. WOW. Aug 20 at 16:59
• Great to hear that! Good luck for your project! Aug 20 at 17:12

May be you can deal with the implicit equation for a "conical, elliptic cylinder" x^2 + aa y^2 == Abs[(z bb + cc)]. Gained a factor of three in speed with my PC.

SetDirectory[NotebookDirectory[]];
points = Import["points.csv"];

gp = Graphics3D[Point@points, Axes -> True];

points2 = points /. {r_, s_, t_} -> {r, s, t, 0};

(ff = FindFit[points2 // Rationalize[#, 0] &,
x^2 + aa y^2 -
Abs[(z bb + cc)], {{aa, 11/10}, {bb, 2/100}, {cc, 4/10}}, {x, y,
z}]) // AbsoluteTiming

(*   {0.0155762, {aa -> 1.04171, bb -> 0.00262419, cc -> 0.25958}}   *)

cp = ContourPlot3D[
0 == ((x^2 + aa y^2 - Abs[(z bb + cc)]) /. ff), {x, -1, 1}, {y, -1,
1}, {z, -1, 1}, Mesh -> None];

Show[cp, gp] • Yeah, that's where I started. Your code will turn into my code if you add xy-translation and three rotations :) Aug 20 at 8:35