# Ellipsoid problem

This topic is 3867 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

I am trying to work out a way to draw a line along the surface of the ellipsoid from two known points on the surface. I am planning to divide the line between two known points on the surface into sections and then project the newly divided points inside the ellipsoid (due to the curvature) to the surface If you take the equation of an ellipsoid defined by x^2 y^2 z^2 --- + --- + --- = 1 a^2 b^2 c^2 If you have a point defined as xyz below the surface, I know from the point you can determine the direction from the origin, but how do you determine the remaining distance to the surface from the point below the surface? I know that the following equation will give me the normal tangent of the ellipsoid surface, but I dont think this is what i need since I am assuming xyz needs to be on the surface to compute the tangent normal vector? 2x^2 2y^2 2z^2 ---- + ---- + ---- = Normal a^2 b^2 c^2 Could someone offer advice Regards James

##### Share on other sites
The hull is described by a implicit function
f(x,y,z) = 0
where in the case of the ellipsoid
f(x,y,z) := x2 / a2 + y2 / b2 + z2 / c2 - 1

A line passing through 0 and also through an arbitrary point p is descibed by a vector valued function
l(k) := 0 + k * ( p - 0 )
where all those nasty 0s can be omitted in computation but are shown here for completeness, and to retain the sematics of point vs. direction vectors.

Now, l can be expressed by 3 scalar functions as well
x(k), y(k), z(k)
and these can be set into f(x,y,z), leading to a quadratic function with just k as unknown variable. In this case, solving for k yields always in 2 results
k1, k2
as long as p isn't 0 itself, but the negative k (let's asume it is k2) is to be dropped as long as you don't want the hit on the opposite side, too.

Setting in the remaining k1 into l(k) gives the point of interest
l(k1)
where the line hits the surface. The distance to the surface (measured from the origin) is then
d := | l(k1) - 0 | = k1 * | p - 0 |

Notice please that the direction ( p - 0 ) is in general not identical to the direction of the normal on the surface when looking at the hit point!

[Edited by - haegarr on March 19, 2008 2:22:15 PM]

##### Share on other sites
Thanks for your reply, am going to knock up a quick opengl test harness to try out the theory and see how I get on, tend to find it eariser to visualize and understand things on screen than on paper.

I suspect that there could be another approach using a parameterized curve equation based on time? I am pretty new to 3d modelling and want to learn and understand various techniques that can be employed. I have given myself a little project to try and create a texture mapped google earth type application! ;o) lol

I am assuming that that k1 and k2 are Quadric coefficients?

Cheers

[Edited by - osso on March 19, 2008 6:53:57 PM]

##### Share on other sites
Times like these I wish I paid more attention in maths alevels 10 years ago!!

Where is a good place to read up on quadratic equations?

this is how far i've got up to now.

given the line l(k) := 0 + k * ( p - 0 )

x = x0+xd*k
y = y0+yd*k
z = z0+zd*k

subsitute these back into the ellipsoid equation

(x0+xd*k)^2/a^2 + (y0+yd*k)^2/b^2 + (x0+xd*k)^2/c^2 = 1

this is where i start to struggle with quadratic equations. I assume that I need to rearrange and simplify the equation for k so it looks somewhat like this?

A*k^2/a2 + B*k/b2 + C/c2 - 1 = 0

Regards

##### Share on other sites
You're on the right way. Notice please that 0 means a vector valued zero. So, from a pure mathematical point of view, the formula
l(k) := 0 + k * ( p - 0 )
could be written simplified as
l(k) := k * p
I've chosen to use the full formula above because it retains the semantics
point + scaled direction => point
of vectors. In 3D graphics, sometimes the lack of knowledge of the difference between point vectors and direction vectors causes some confusion (several threads here on GDnet attest this fact).

The line formula chosen by you is more general, since it allows to work with lines passing through any "origin" [x0 y0 z0] instead of just through zero. We can follow this approach and simplify the solution later on. But you can set all of the x0, y0, z0 to zero from the beginning if you want to make your life easier. I definitely propose you to verify my math at all, because it is so easy to write down a bug in such a post :(

The substitution is then (you've found it correctly)
( x0 + k * px )2 / a2 + ( y0 + k * py )2 / b2 + ( z0 + k * pz )2 / c2 - 1 = 0

With the definition
m := a2 * b2 * c2
we can multiply all terms with m, yielding in
( x0 + k * px )2 * ( b * c )2 + ( y0 + k * py )2 * ( a * c )2 + ( z0 + k * pz )2 * ( a * b )2 - m = 0

Looking at the parts containing k, we have to resolve the brackets. E.g. for the "x" term
( x0 + k * px )2 * ( b * c )2 = k2 * px2 * ( b * c )2 + 2 * k * px * x0 * ( b * c )2 + x02 * ( b * c )2
Doing so for the "y" and "z" terms also, and collecting all resulting terms with k2, k1, and k0==1, resp., we get something of the form
r * k2 + s * k + t = 0
or, after a simple divison
k2 + s/r * k + t/r = 0
where r, s, and t are relative complex terms, at least if you haven't used the 0 origin.

Now something comes into play that is called the "p,q formula" here in Germany. I don't know how it is called in your country. However, the formula says that the general solution to the equation above is
k1 := -s/(2r) + sqrt( ( s/(2r) )2 - t/r )
k2 := -s/(2r) - sqrt( ( s/(2r) )2 - t/r )

It may happen that the term under the square root is positive, in which case 2 (real) solutions exist, i.e. k1!=k2. It may happen that the term under the square root is 0, in which case 1 solution exists, i.e. k1==k2, and it may happen that the term under the square root is negative, in which no (real) solution exist. To interpret the results: 2 (real) solutions mean that the line hits the ellipsoid in 2 distinct places. 1 solution means that the line tangents the ellipsoid only, and 0 (real) solutions mean that the line misses the ellipsoid.

In your special case, where the "origin" of the line is 0, I expect that
k1 == -k2
will be the results. A negative k means that the line hits the ellipsoid when "traveling backwards" on the line, and hence against the direction in which p lies w.r.t. the global origin, so I had said that the negative k will be no solution you are interested in.

EDIT: Made the stuff a bit easier to understand by dropping some details.

[Edited by - haegarr on March 20, 2008 5:43:09 AM]

##### Share on other sites
Haegarr thanks again! I managed to come up with a general quadratic solution for k after my post and came up with the following on paper

k = (-b + sqr(pow(b,2) - 4ac))/2a
k = (-b - sqr(pow(b,2) - 4ac))/2a

which appears to be rather similar to your solution, however your post does help to clarify a number of points so thanks again for your time. I will try and convert to opengl code and plot the intersection points on the ellipsoid.

I suspect this is not the most elegant way of projecting a line along the surface given a point on the surface, a direction, and a distance you with to travel along the surface of the ellipsoid. I have heard of parameterized curve equations based on time but at this stage i just dont know enough about them.

Kind regards

##### Share on other sites
Quote:
 Original post by ossoI suspect this is not the most elegant way of projecting a line along the surface given a point on the surface, a direction, and a distance you with to travel along the surface of the ellipsoid. I have heard of parameterized curve equations based on time but at this stage i just dont know enough about them.
For a general hull, it it a good solution. You can try to compute the curvature on the surface to select reasonable sample point distances, so that errors can be constrained to an upper limit.

However, just for furthur thinking, especially for an ellipsoid there is also another way: You can find a parametric representaion of the curve on a sphere, and them scale the curve (i.e. transform its base points) accordingly to the space transformation needed to map the sphere to the given ellipsoid. If you think about those scaling transformation, it is just the inverse you already do. With
x2 / a2 + ... == ( x / a )2 + ...
where a is the radius in principal x direction, you scale the space with the ellipsoid by 1/a along the principal x axis. Doing so for all 3 axes means to scale the space so that the ellipsoid is made a (unit) sphere. This, of course, works only if the ellipsoid is aligned with the said axes, but that can ever be yielded by a rotation.

##### Share on other sites
One thing that might help to solidify things in your mind is to note that the ellipsoid and the hyperboloids form a triply orthogonal system.

If you generate a hyperboloid and ellipsoid, you will see that the surface normal of the ellipsoid where it intersects with the hyperboloid is also a tangent vector of the hyperboloid.

The first image on this page shows the intersection of an ellipsoid and two hyperboloids (one of one sheet, one of two sheets):

http://www.math.uni-bonn.de/people/weber/research/triple/index.html

The book Modern Differential Geometry of Curves and Surfaces with Mathematica by Alfred Gray also discusses this topic. It's pretty cool stuff.

If you search on Google for "triply orthogonal differential mathematica", you should get a preview of the book. The spherical equivalent of the triply orthogonal system is on page 606 (with another example of the ellipsoidal version on page 615).

I have the 2nd edition of this book and there was a whole chapter dedicated to the topic. Apparently this 3rd edition only spends part of a chapter on it. Strange.

Anyway, the day I saw these images, I immediately fell in love with the ellipsoid/hyperboloids.

[Edited by - taby on March 20, 2008 5:10:09 PM]

##### Share on other sites
Quote:
 Original post by haegarrHowever, just for furthur thinking, especially for an ellipsoid there is also another way: You can find a parametric representaion of the curve on a sphere, and them scale the curve (i.e. transform its base points) accordingly to the space transformation needed to map the sphere to the given ellipsoid.

I'm with you on this.

To start, I'm going to offer some more compact notation. This roughly summarizes what you said, but I think you'll agree that it's pretty this way. :-)

Matrix notation

Any ellipsoid can be expressed by the implicit function,

x' P x = 1

where x is a vector (in 3-space, in our case, but this works in any number of dimensions), and P is a positive definite matrix. I use x' to mean "x, transposed."

The eigenvectors of P give you the principal (major and minor) axes, and the eigenvalues give you the squares of the reciprocals of the lengths of the ellipse in these directions. If you diagonalize P, this is,

x' V Y V' x = 1

where the columns of V are (orthonormal) vectors along the principal axes, and Y is the diagonal matrix of eigenvalues -- or, of the squares of the reciprocals of the axis lengths.

Change of coordinates

If we define z = V' x, then this becomes,

z' Y z = 1

This is just the equation for the same ellipsoid in another coordinate system. Really, we have just rotated the ellipsoid so that its principal axes line up with the coordinate axes. The transformation from the original coordinates to the new coordinates is given by V'.

We can make this even better. Let D be a diagonal matrix whose elements are the square roots of the elements of Y. That is, D contains the reciprocals of the axis lengths. Then,

Y = D' D

so our ellipse equation becomes,

z' D' D z = 1

and defining,

w = D z

then our ellipse equation reduces to,

w' w = 1

Note that w' w is just the sum-of-squares of the elements of w; i.e., it is the squared 2-norm of w. So this equation says, "the length of w is a constant: 1." In other words, it describes a sphere.

So, what have we done? By transforming x to w, we have changed our ellipse into a sphere. To summarize, that transformation was,

w = D V' x

and the inverse is,

x = V D^(-1) w

where D^(-1) just contains the reciprocals of the elements of D -- or the axis lengths! For an even more compact notation, let's just define,

T = D V' <==> T^(-1) = V D^(-1)

so,

w = T x

Can't get simpler than that!

Drawing the arc

Suppose we want to get from x1 to x2 on the surface of the ellipsoid. We will equivalently go from w1 to w2 on the surface of the unit sphere, where w1 = T x1, w2 = T x2.

There are two ways to do this. I'll start with what I think is the easier-to-understand way. (But the second is probably easier to implement).

Subdivision

It's very simple. Given two vectors w1 and w2 on the unit sphere, we can compute their midpoint and then normalize it, to get another point on the sphere halfway between w1 and w2:

wm = normalize{ (w1 [plus] w2)/2 }

By "halfway between," I mean that the angle from w1 to wm is the same as the angle from wm to w2, so the arclength along the sphere from w1 to wm is the same as that from wm to w2.
(Side note: Can anybody tell me why my plus signs don't show up? That's why I used text above.)

If we do this recursively with a recursion depth n, we can obtain 2^n[plus]1 evenly-spaced points along the arc you want on the sphere. Transform them to x-space using T^(-1), and voila, you have a bunch of points on your ellipse along the arc you want.

Since this method doesn't require trig functions -- just some square roots -- I'd expect it to be reasonably fast.

[The one time this method breaks down is when w1 = -w2 (in which case the midpoint is the origin, and if we try to normalize we divide by zero). Here, any arc from w1 to w2 is the shortest path along the surface, so we can pick one arbitrarily. To do this, we need a vector -- any vector -- perpendicular to w1 and w2 which we can use as the midpoint for the first iteration of the algorithm; after this, we can subdivide normally. An easy way to do this would be to generate a vector at random (with elements uniformly distributed between -1 and 1, perhaps), subtract off its projection onto w1, and then normalize it.]

Parameterization

There are a lot of ways to do this. Here's one (and probably the best one, at that): SLERP! Wikipedia's article on the subject gives a very nice, simple formula for doing this: http://en.wikipedia.org/wiki/Slerp

After you SLERP in w-space, again, just transform by T^(-1) to get back to x-space.

The thing I like about this solution (which it has in common with the subdivision solution) is that it generalizes to an arbitrary number of dimensions.

There you go. Hope you enjoy!

1. 1
Rutin
45
2. 2
3. 3
4. 4
5. 5
JoeJ
19

• 13
• 10
• 12
• 10
• 13
• ### Forum Statistics

• Total Topics
632999
• Total Posts
3009814
• ### Who's Online (See full list)

There are no registered users currently online

×