# Constraining a least squares minimization to fit a single root?

I have an algorithm that segments depth images using surface fitting. At the moment the algothim uses least squares polymonial fitting, but polynomials are not powerful enough to fit the shapes that are in these images. I replaced the explicit polynomial $z = f(x,y)$, with the implicit fitting problem $f(x,y,z) = 0$. Where $x$ and $y$ are the pixels location and $z$ is the value of the pixel and $f$ is a polynomial. A classic problem with many nice linear solutions. This did make the fitting much more powerful, but this left me with a problem I have been unable to solve for a while now. Once I had the least squares solution to the implicit poly, I had to solve for z, not easy at all, but I could search for the minimum (there are only 256 possible pixel values to search through after all). THE problem is that there was more than one solution for z! Where a pixel can only have one value. That is, once the implicit poly was solved for z there were several roots.

So my question is;

How do I formulate a least squares minimization problem for a single root? What contraints must I add?

This might not be completely clear so here is an example: You have a set of noisy data points $\{x,y\}$ that form a semi-circle about the origin in the positive y only. I want to fit $f(x,y) = 0$ to the data points, or to be pricise, a circle. $f == x^2 + y^2 + c = 0$.

The least squares minimization problem is a nice simple linear $\min_c \sum_{i=0}^n (x_i^2 + y_i^2 + c)^2$ but I am only interested minimizing the datapoints distance to the positive half of the circle. Solving for $y$ and taking the positive root gives us $y_i = \sqrt{-c-x_i^2}$. Giving us the acutal minimization problem as $\min_c \sum_{i=0}^n (y_i - \sqrt{-c-x_i^2})^2$ NOT a nice linear problem at all. These two minimisation problems will give differnt results since the negative half of the circle should not try to fit itself to any data points. How can i constrain the linear problem so that it is an equivelent to minimizing to a single root?

Or in general constrin $\min_{ijk} (\sum_{i,j,k}a_{ijk}x^iy^jz^k)^2$ to a minimization to a single root in $z$?

I hope this makes sence and sorry for the length but this has been annoying me for weeks.

-

Asuming that your surfaces are not too wild, it seems logical to try to combine the implicit approximation with the explicit one, i.e., to solve 2 least square problems simultaneously: one for the implicit function and one for the explicit one. The explicit solution will give you a crude approximation for the root and then you just choose the root of the implicit function that is closest to this crude approximation as the "true" position of your pixel. It certainly works with your half-circle example but whether it'll work in your real case depends on the actual surfaces you are dealing with.

If you could give us some idea of what they look like, what precision you are getting with the explicit approximation, what precision you are getting with the implicit one, and how close the "parasitic" roots can come to the "true" ones, we might be able to say more.

OK, let's try one more idea. Minimize the sum of $f(x,y,z)^2+A(\partial_z f(x,y,z)-1)^2$ with some positive $A$ (you'll have to play with its choice to see what works best). The advantage is two-fold: first, now you have the true root separated (the $z$ derivative is not small) and second, you have two approximate equations for it ($f=0$ and $\partial_z f=1$), which should give you extra advantage, the hope being that the parasitic roots of $f$ won't be able to match the derivative too. Also, keep the power in $z$ well below that in $x$ and $y$. Note that you are still better off than with $z-f(x,y)$ because this explicit formula satisfies all the extra conditions we are trying to impose automatically.

And yes, it'll help to see the data, though, if possible, I'd prefer the one-dimensional case (the $x$ slice of your real data should have all the same problems already). Just post something reasonable (like $S(1),\dots, S(50)$) that we want to approximate by $f(x,z)=0$ with $f$ of some reasonable degree with which you currently have a problem. I'll try to play with it a bit when I have free time.

-
I like this idea, but the "parasitic" roots can get very close to the root i am interested in (esp at high polynomial degrees) and this solution does not stop these roots from trying to fit themselfs to the data. I will post some images of what I am trying to do here if that will help, but I am glad to see that you understand the problem. – Tategami Mar 11 '10 at 10:09
I hope I do (I edited my post to add one more suggestion, by the way). But it is hard to see what exactly you are against with without seeing real data... – fedja Mar 11 '10 at 13:45
I must admit that I dont quite understand your second solution. You are trying to contratin the z gradient to a constant at the roots? How does this constrint apply to one root and not the other? I should have some examples ready for you soon. Thanks – Tategami Mar 12 '10 at 15:51
Not at all the roots! You minimize the sum over the actual set of (x,y,z) values, which is one particular point for each x and y, not all roots. So, at the point that you really want both the value and the gradient are expected to be close to the prescribed values but everywhere else they are completely unrestricted. – fedja Mar 12 '10 at 18:39
Does my example below help you see what I am trying to do and make my question clearer? The other side of the problem is that I also dont know which root has been fitted and so I dont know which minimum to search for. To see if I understand your second suggestion, this is how I understand it: You are adding a type of regularisation term, by forcing the gradients in the z to be the same (or at least close) we are making the surface smoother. Further to this, one might expect different roots to have different gradients in z. I hope to have time next weekend to try this solution. Thanks. – Tategami Mar 22 '10 at 12:53

I think that it is better to use rational approximation for the fitting z = P(x,y)/Q(x,y), where P and Q are polynomials. This type of fitting is more flexible in general than polynomial fitting and doesn't have the multiple solution problem. There are many algorithms to establish this type of fitting, one possibility is the Nelder-Mead simplex algorithm.

-
I have thought about rational approximation but I don't think it will do what I want. What kind of shapes can I fit with this method? How well will it do with the circle example that I give above? Thanks – Tategami Mar 11 '10 at 20:49
In general, rational functions can fit a wider range of shapes than polynomials, see for example the following Wikipedia page, where the two types of fittings are compared en.wikipedia.org/wiki/Polynomial_and_rational_function_modeling – David Bar Moshe Mar 13 '10 at 15:38

I finally have an example ready...

So we have an image of a sphere, and we take data points on it.

Since this is a sphere it is reasonable to assume that we can fit a sphere to this part of the image. Here are the reults:

Blue is low error, we can see that most points fit well, the other are the `parasitic' roots, interpolations look as follows:

The correct root fits the image well, the other root is just getting in the way (at the moment) but here is the fitted function:

Now this is a nice result since the 'edge' of the fitted sphere corresponds to the edge of the sphere in the image. Something I imagine the suggested rational fitting can't do any why I WANT several roots in the solution, but I only want to fit one to the data!!

So whats the problem, well, this is a nice result but a 3rd order polynomial might fit it better (the data is technically not a simple sphere), or i might have other data that would need even higher order fitting (toruses for exmaple are 4th order). So I fit a 3rd order poly and here is the resulting function:

Not nice, the roots are all over each other, fighting against each other to fit the data.

Now I want the re-iterate my question. One can always solve for z and fit a single root to the data using non-linear regression, but my feeling is that we should be able to do the same using constrained linear fitting. And HOW is my question! Here is an example that constraining the fit does at least help. Here the $z$ and $z^2$ parameters are fixed to $1$: Fitting the same 3rd order poly to the same data:

A much friendlier 3rd order poly, but this restriction is too much it turns out, and i have lost the nice results.

So the question in terms of my original example:

How can I linearise the following minimization problem without reintroducing multiple roots?

$\min_c \sum_{i=0}^n (y_i - \sqrt{-c-x_i^2})^2$

and then we can see if this can be done in general to all implicit polys....

Thanks.

-