0
101 Nov 05, 2007 at 17:43

To determine the closest point on an ellipsoid’s surface to any given other point (x0, y0, z0), I have found a method of parameterizing the ellipsoid equation by a parameter t, so that:

f(t) = a² * x0² / (a² + 2t)² + b² * y0² / (b² + 2t)² + c² * z0² / (c² + 2t)² - 1

Which you will find is simply the standard equation, where x, y and z have been substituted by:

x(t) = x0 * a² / (a² + 2t)
y(t) = y0 * b² / (b² + 2t)
z(t) = z0 * c² / (c² + 2t)

Since the 1 has been pulled to the other side, the root of f(t) will be the value for t at which the components x(t), y(t) and z(t) of a new point will be such that the new point is on the ellipsoid surface. (Since if the equation is 1, this is the case.) This point is the one we are searching for.

To find this root, a simple Newton-Raphson approximation is applied to the function. I’ve tested this out, and it works just fine. However, I’m not sure why x(t), y(t) and z(t) are defined the way they are.

It is said in the source of this approach that it is done this way because “the derivative of the ellipsoid equation at the nearest point on the ellipsoid to the skin point represents a vector between the skin point and its nearest point”. This is one thing I don’t grasp. The other is that the three functions are called components of a “parametric line function”, though when I input different values for t and visualized this, the resulting path in 3d space was clearly a curve, not a straight line.

The author then just describes the approximation as moving step by step on the line/curve created by the three functions towards the t value that gives a point almost on the ellipsoid surface. Another important property of the functions he mentions is that for t = 0, the resulting point on the line/curve is exactly the point we want to find the nearest point for (as one can easily see).

So, all in all, I’d like someone to tell me how the author got to these three equations.

#### 23 Replies

Please log in or register to post a reply.

0
165 Nov 05, 2007 at 19:04

I don’t think those equations are right… I just tried them now in the 2D case, and they didn’t give the right answer at all. (I used an ellipse with a = 2 and b = 1, and (x0, y0) = (2, 1.5), and plotted the parametric curve according to the equation you posted. It goes through a point on the ellipse, but definitely not point on the ellipse nearest to (2, 1.5).)

But, assuming that some typo was made, I believe I can give a partial explanation:

Imagine taking the point (x0, y0, z0) and drawing a small sphere around it, then gradually expanding that sphere till it just touches the ellipsoid. The point where it touches the ellipsoid is the place you’re looking for, and since the sphere and ellipsoid are just touching at that point, they’re tangent there. This means the ellipsoid’s normal at that point is in the opposite direction to the sphere’s normal - and since the sphere’s normal points directly away from its center, it follows the ellipsoid’s normal points directly towards the sphere center, aka (x0, y0, z0).

So we know the location we want is the point on the ellipsoid where its normal points toward (x0, y0, z0). Now, if you define an ellipsoid as a scalar field f(x,y,z) where f = 0 for points on the surface, f < 0 inside and f > 0 outside, then the normal of an ellipsoid is given by the gradient of f.

Now, here is where I’m not sure what your author is doing anymore - he seems to be trying to construct a 1D curve that (somehow) will go through the right point on the ellipsoid. At this point it makes more sense to me to parameterize the ellipsoid as a function P(u,v) mapping u,v to some point on the ellipsoid surface. Then, what we want to find are the u,v for which

(x0, y0, z0) - P(u,v) = k * gradientF(P(u,v))

for some k, i.e. the vector from the point on the surface to (x0, y0, z0) is in the same direction as the gradient at that point. This can be re-expressed as a function that normalizes both vectors and takes their dot product:

G(u,v) = normalize((x0, y0, z0) - P(u,v)) dot normalize(gradientF(P(u,v)))

Now, this function will have a maximum exactly at the u,v coordinates where you’ll find the point you’re looking for. You have to apply a two-dimensional maximization now, but that’s not so bad as the function is very well behaved, there are no local maxima. You could even probably solve it analytically, though it would be messy (and possibly the numerical solution would perform better).

0
101 Nov 05, 2007 at 22:33

Well, they don’t seem to be right if they should be doing what the author claims them to do, although I’ve read a few more of his papers and almost all of them had errors of some kind in their equations.

Also I’m not sure the procedure really gives you the closest point of the surface, but it certainly gives some point on the surface, so when I quickly tested this yesterday, I thought it to be working all right. That might have been wrong, I don’t know.

Concerning the equation f(t), this is nothing different than 0 = x²/a² + y²/b² + z²/c² - 1, with x, y and z replaced by the three functions x(t), y(t) and z(t). So if you had a typo in there, you should be able to find it.

And now to your explanation: I understand what you’re saying with the normals of the ellipsoid and the imaginary sphere (and expanding your sphere would be nothing else than moving along the line described by the author). But he talks about the derivative of f(t) at the nearest point representing a vector between that point and the reference point. But f(t) returns a scalar value, so how can this be any kind of vector or even any kind of normal?

I think there is just some consideration behind these component functions that I just don’t get, though I think it might be pretty simple.

Well, to be honest, the second part of your post was pretty cryptic to me (because I’ve never had any vector analysis in school), though I understand the part with the gradient, which makes sense to me (even though how to calculate it is a whole different story).

I’ll work some more on this tomorrow, but now I unfortunately haven’t got the time. Even though your solution, as little as I understand of it, makes a lot more sense to me right now, I’d still like to understand what the author of this paper was thinking. After all, he probably had a point in doing it his way.

0
165 Nov 06, 2007 at 00:27

@Sniffman

But he talks about the derivative of f(t) at the nearest point representing a vector between that point and the reference point. But f(t) returns a scalar value, so how can this be any kind of vector or even any kind of normal?

I think he must mean f(x,y,z), the ellipsoid’s implicit function. This is a function from vectors to scalars, and the derivative of such a function is a vector, called the gradient, which is defined as the vector of the partial derivatives of the function with respect to each axis. So it’s (df/dx, df/dy, df/dz). This vector is normal to the ellipsoid.

Of course, when you substitute x, y, z with x(t), y(t), z(t) you get f(t) which is a function from scalars to scalars, which has only a regular scalar derivative (which is the dot product of the gradient of f(x,y,z) with the tangent vector to the curve (x(t), y(t), z(t)), by the chain rule).

I’d still like to understand what the author of this paper was thinking. After all, he probably had a point in doing it his way.

I’d like to understand it to. If there really is a way to do this using 1-dimensional root finding rather than 2-dimensional maximization as I described, then it would probably be better. Can you give a link to the original paper? Maybe I can make some more sense out of it.

0
101 Nov 06, 2007 at 18:06

A quick glimpse using Google tells me you should be able to download the paper here. The topic of interest is chapter 3.7. I’m warning you, though: There’s not much of an explanation in there.

Just as a test whether the two of us can reproduce results with the methods used so far. What’s the point you get for x0=1.5, y=1.5, z0=0.5, a=2, b=1, c=1? For these values, the root of f(t) is at t=0.490372, which gives me approximately the point (1.2, 0.7576, 0.2525). Can you confirm this?

Assuming you now read the paper, I just noticed that I might have gotten wrong what the author calls a line. Seems he doesn’t mean a line of points (x(t), y(t), z(t)), but rather the line dt. But that is nonsense, since the last two equations are simply the Newton-Raphson steps of dividing the function through its derivative and subtracting this from the current value. But it could also mean taking steps OF dt along the line. It’s these unclear expressions that make this so hard to understand.

Well, maybe now you have an idea why x, y and z are chosen that way.

0
165 Nov 06, 2007 at 19:42

@Sniffman

Just as a test whether the two of us can reproduce results with the methods used so far. What’s the point you get for x0=1.5, y=1.5, z0=0.5, a=2, b=1, c=1? For these values, the root of f(t) is at t=0.490372, which gives me approximately the point (1.2, 0.7576, 0.2525). Can you confirm this?

Yes, those are the numbers I get, too. I graphed the ellipsoid and the point, and it does look like it’s the correct nearest point (maybe the equations should vary according to dimension, and that’s why it didn’t give the right answer for the 2D case I tested?)

Well, maybe now you have an idea why x, y and z are chosen that way.

You were right, there’s really no explanation in the paper. I’ll keep thinking about it, but for now I’m stumped.

0
101 Nov 09, 2007 at 21:45

I played around with some of your considerations. If you take a sphere around the point again and expand this so far as to be tangent with the ellipsoid, then the normal of the sphere will be the opposite of the normal of the ellipsoid. The normal of the ellipsoid itself is the gradient of its implicit function.

In other words:

n_sphere = normalize((xt, yt, zt) - (x0, y0, z0))
n_ellipsoid = gradient(x²/a² + y²/b² + z²/c²)
n_ellipsoid = (2 x/a², 2 y/b², 2 z/c²)

then there is:
n_ellipsoid = -n_sphere

which is:
(2 xt/a², 2 yt/b², 2 zt/c²) = normalize((x0, y0, z0) - (xt, yt, zt))

This is correct, as I have tested this for the 2d case and it works. Calculating the gradient and mapping its normalized opposite at the right position on the ellipse boundary will give me what appears to be a normal. The normal vector is also correct for a sphere around (x0, y0, z0), as it will pass through the center point if expanded.

That’s also about what you wrote some posts earlier, just that you introduced a certain k. Was this a vector or a scalar, so was your multiplication a scalar one or a dot product? And what’s with the additional function f(u, v)? Why not just solve this after (xt, yt, zt) and be done?

0
165 Nov 10, 2007 at 01:08

It was a scalar. All I was saying is you want the point on the ellipsoid such that the normal there points toward (x0, y0, z0). If A, B are vectors, A = k*B can only be true if A and B are parallel, otherwise there’s no scalar that can be multiplied by one to get the other.

I’m not sure what you mean by “solve for (xt, yt, zt).” The problem is that since you don’t know what the point on the ellipsoid surface is, you don’t know what the normal should be. Instead you want whichever point happens to have a normal that goes in the right direction. So, I suggested parameterizing the ellipsoid, which assigns each point on its surface a (u,v) coordinate, then constructing a function G(u,v) that will attain its maximum value exactly at the (u,v) where the ellipsoid’s normal is in the same direction as the vector from the ellipsoid point to (x0, y0, z0). (This works by normalizing the two vectors and taking their dot product. We know the dot of two normalized vectors attains its greatest value, 1, when the vectors point in the same direction.)

Actually, now that I think about the problem again, I have another idea. If you think about the sphere again, the ellipsoid’s function produces a positive value on each point of the sphere. (Recall the ellipsoid is defined as a scalar field f(x,y,z), which is zero at points on the ellipsoid surface, < 0 for points inside, and > 0 for points outside.) Now, there will be some point on the sphere at which the ellipsoid function takes a minimum value. If you can solve for the location of this point, you could parameterize it by the radius of the sphere, and thereby obtain a parametric curve that starts at (x0, y0, z0) and goes toward the ellipsoid, intersecting it at the nearest point. It probably wouldn’t be a straight line, but if you imagine the sphere growing till it just touches the ellipsoid, then at the point where they touch, the ellipsoid function would have value 0, while at other points on the sphere it would have value > 0, so the point where they touch is definitely a point on the parametric curve. This might be the function that your author used, but I haven’t derived the equations for it, so I don’t know for sure.

0
101 Nov 11, 2007 at 17:38

I see what you mean. That might be a way, but, as you already noted, it seams far away from what is done in the paper.

What I meant was simply to re-express the equation so that it reads (xt, yt, zt) = something. Now, I’m not very experienced with handling vector equations, but I know that would be a strategy when dealing with scalars. This way, I thought one might be able to express (xt, yt, zt) to be somehow depending on (x0, y0, z0) and the ellipsoid dimensions.

As for your other consideration: This doesn’t sound bad. But how to find the point on the sphere for any arbitrary radius that will produce the minimum ellipsoid value of all possible points on the sphere with that radius?

0
165 Nov 11, 2007 at 19:05

Well the obvious way is to parameterize the sphere and use the usual technique for finding the minimum of a function (take the partial derivatives of the ellipsoid function with respect to the sphere parameters, set equal to zero and solve). That’s messy, but if you work it all the way through you might be able to get something simple at the end. There may well be something I’m missing and a simpler way to do this. I don’t have time to work out all the details right now, though.

0
101 Nov 11, 2007 at 19:59

I didn’t read through the post and maybe this was suggested already, but because you get ellipsoid by applying linear transformation on a sphere, you can transfer the point with inverse ellipsoid transform, compute the closest point of the sphere and transfer the point back.

0
165 Nov 12, 2007 at 00:31

Unfortunately, that doesn’t work since the transform does nonuniform scaling and therefore distorts distances.

0
101 Nov 12, 2007 at 06:25

You don’t compute the distance after the inverse transform, but the closest point on the sphere and transform that point back. Then compute the distance between the two points in the “non-distorted” space. Sorry for not being explicit enough in my post.

0
165 Nov 12, 2007 at 07:33

No, I still don’t think it works. The closest point in the distorted space (where the ellipsoid is a sphere) doesn’t transform to the closest point in the non-distorted space.

Think about it: the set of points of distance r from a point is a sphere of radius r in the non-distorted space, but in the distorted space that would be transformed into an ellipsoid, so we’d be back to the same problem again.

0
101 Nov 12, 2007 at 19:30

Okay, I think I’m starting to understand what you want, but it’s an awfully slow process. Thing is, I’ve never done such things as partial derivatives and such stuff. So I’m obviously one step behind you. I have a rough idea what you are trying to do, but it would be nice if you could just list the steps you want to take.

Right now, I’ve taken you as follows: You want to create the partial derivatives of the ellipsoid equation (which would be 2x/a, 2y/b and 2z/c) and parameterize these equations by exchanging x, y and z in these equations with something including the radius of the sphere around the point. Then, you want to find the minima of these functions. This is where I’m having my problems. What do you solve this after? And how will this help us? Or rather: How does this translate ultimately into functions for x, y and z depending on some parameter t, as in the original paper?

0
165 Nov 13, 2007 at 03:59

Okay, forget about my post #10 above, I just remembered something else: Lagrange multipliers. This can make our lives a lot easier - it’s a mathematical technique for finding the minimum and maximum of a function given a constraint, which is exactly what we want to do.

Suppose the ellipsoid function is f(x,y,z) = x\^2 / a\^2 + y\^2 / b\^2 + c\^2 / z\^2. (Actually, there should be a - 1 at the end there, but I’ve left it off to make things a bit simpler since it doesn’t affect where the minimum and maximum are.) For a given radius r, the sphere function is g(x,y,z) = (x - x0)\^2 + (y - y0)\^2 + (z - z0)\^2 - r\^2. Now we want to find the point (x,y,z) that minimizes f(x,y,z), subject to the constraint g(x,y,z) = 0. The theory of Lagrange multipliers tells us how to do this. First we construct the function:

L(x,y,z, lambda) = f(x,y,z) + lambda * g(x,y,z).

Here, lambda is an additional scalar, so the function is now four-dimensional. Then we need to find the point(s) where the gradient of this function is zero - in other words, where all of its partial derivatives are zero. So, we have to take partial derivatives of L and set them all equal to zero, producing a system of four equations. Most likely we can solve these equations analytically. There should be two solutions: one for the maximum and one for the minimum; we throw away the maximum one since we don’t care about it. The solution will express x, y, z, and lambda as functions of x0, y0, z0, a, b, c, and r. We can throw away the equation for lambda as it’s not important anymore. The other three equations will tell us the (x,y,z) point that minimizes f on the sphere of radius r about (x0, y0, z0). Now if we vary r, we will produce a curve in space that (as discussed before) is guaranteed to go through the desired point on the ellipsoid. So all we have to do is substitute the curve equations x(r), y(r), z(r) into the ellipsoid function f(x,y,z) to get f(r), and find the root of that function using Newton’s method.

0
101 Nov 13, 2007 at 05:07

@Reedbeta

No, I still don’t think it works.

Yeah, you are probably right. Just some tongue-in-cheek thinking here: Is the normal of the sphere where it intersects the line from origo to the point the same as the normal of the closest point in non-distorted space ellipsoid? If that’s the case, you should be able to find the point with transforms.

0
101 Nov 13, 2007 at 17:37

Okay, I started doing this. My function is:

L(x, y, z, lambda) = x²/a² + y²/b² + z²/c² + lambda ((x - x0)² + (y - y0)² + (z - z0)² - r²)

The partial derivatives are easy to get to, but I’, having a question regarding the derivative for lambda. This is of course exactly the sphere equation. But lambda itself doesn’t appear in this equation, so it’s impossible to solve after lambda (though I believe you said that we don’t have to). Also, if that is correct, I do need to substitute lambda in the three equations for x, y and z by the partial derivative I get for it, right? But if that is the case, it seems rather unlikely I’ll be able to get most of the variables and constants out again. For instance, the equation in the original paper for x was only dealing with t, x0 and a, and the other coordinate components behaved accordingly. I don’t see this happening if what I think I have to do is correct.

Is there some apparent strategy for solving this system of equations?

0
165 Nov 13, 2007 at 20:42

Lambda won’t appear in the equation gotten from taking the derivative for lambda, but it will appear in the other three equations.

Remember, you’re not substituting lambda by the partial derivative. You take each the partial derivatives and set it equal to zero, producing a system of equations that has complicated stuff involving x, y, z and lambda on one side and all zeros on the other side.

It’s going to be nonlinear, so no, there’s no specific strategy for solving it, except to solve for lambda in one of the equations (that involves lambda) and substitute it in for lambda in the other equations, thus getting down to three equations that don’t involve lambda. And then repeating the process for another variable. Or you could do a different variable first, and do lambda later. Which order you solve the variables in won’t affect the correctness of the solution, but it may well affect how easy it is to do. You just have to use your intuition and ultimately trial and error.

It’ll be a messy process, but you should be able to get something reasonably simple out at the end. The equation is quite symmetric in x, y, and z so the solution should also be like that.

0
101 Nov 14, 2007 at 21:03

Hmm, so know I have used equation (i) to express lambda by x, and then I have formed equations (ii) and (iii) to express y and z by lambda, but there’s obviously some connection missing here. How do I connect x, y and z with one another?

Right now I have:

lambda = -x/(a² * (x - x0))
y = (b² * lambda * y0)/(b² * lambda + 1)
z = (c² * lambda * z0)/(c² * lambda + 1)

But of course y and z are not connected in any way. Also, I have no idea how to use the former sphere equation to get r into the other equations (because r is of course our most important parameter). Somehow I seem to be missing a point here. I do know that with each new equation used, I should connect two additional variables, but I don’t see how to do this here.

It’s probably pretty easy afterwards, though…

0
165 Nov 15, 2007 at 00:13

Just to check. I got this system (before trying to solve for anything):

0 = 2x/a² + 2 * lambda * (x - x0)
0 = 2y/b² + 2 * lambda * (y - y0)
0 = 2z/c² + 2 * lambda * (z - z0)
0 = (x - x0)² + (y - y0)² + (z - z0)² - r²

As you see, the ‘r’ comes into play in the last equation. Try solving (i), (ii), and (iii) for x, y, and z respectively, and plugging these into the last equation. That should leave you with one equation in one variable, lambda, which can be solved. Then you’ll have lambda in terms of x0, y0, z0, and r. Substitute that into (i), (ii), and (iii) and then solve them again, which will give x, y, z in terms of x0, y0, z0, and r.

0
101 Nov 17, 2007 at 14:37

Yes, these are the equations I’m using, too.

When I re-express the first three equations for x, y and z, I get:

x = (a² lambda x_0)/(a² lambda + 1)
y = (b² lambda y_0)/(b² lambda + 1)
z = (c² lambda z_0)/(c² lambda + 1)

When I input this into the fourth equation, I get:

0 = (x0²)/(a² * lambda + 1)² + (y0²)/(b² * lambda + 1)² + (z0²)/(c² * lambda + 1)² - r²

But now I’m stuck again when trying to solve this after lambda. This will get absurdly long due to the fact that lambda is stuck in theee different binomials in the denominators of three fractions. I don’t even think it is solvable this way.

0
165 Nov 17, 2007 at 18:07

It *might* still be solveable. You can combine all the fractions into one, then since the left side of the equation is zero, you can just set the top of the fraction to zero, and ignore the denominator (though it’s probably a good idea to check what conditions make the denominator zero).

What it looks like is you’ll have a polynomial of degree 6 in lambda. While these aren’t solveable analytically in general, it might have a special form that will allow you to solve it, such as being the cube of a quadratic, or the square of a cubic.

0
101 Nov 17, 2007 at 20:45

Well, of course you are right that you can just leave out the denominator under certain conditions.

But that won’t make the numerator any shorter. And there’s absolutely nothing there to simplify, as every element of the chain of additions we get is unique. So I don’t see any way to get to lambda. Besides, I get a polynomial of degree 4 for lambda. You have to expand the fraction with the other two binomials, respectively. That gives you one lambda squared in each binomial, and if these two are multiplied again you get the fourth degree. How do you get to six?