http://www.geometrictools.com/Documenta ... ipsoid.pdf
this might help
or the retork:
http://www2.imperial.ac.uk/~rn/distance2ellipse.pdf
and/or that:
http://www.devmaster.net/forums/archive ... 10818.html
Closest point to an ellipsoid
Okay, I start to feel sick about ellipsoids right now. I can see nothing else but ellipsoids everywhere...
Thanks for the links, Dorth, the best choice seems to be the second one. It's a bit disappointing though that it talks about x iterations of solving the equations to find that silly point. Finding it will never be as fast and trivial as I first thought.
I'm going to deal with something completely else for a while, just to clear the ellipsoids out of my head. If I won't come up with another idea, I think I'll go for this solution.
But I don't want to "close the topic". Who knows! Somebody might come up with a terrific idea.
Thanks for everyone being involved so far!
Cheers,
PI
Thanks for the links, Dorth, the best choice seems to be the second one. It's a bit disappointing though that it talks about x iterations of solving the equations to find that silly point. Finding it will never be as fast and trivial as I first thought.
I'm going to deal with something completely else for a while, just to clear the ellipsoids out of my head. If I won't come up with another idea, I think I'll go for this solution.
But I don't want to "close the topic". Who knows! Somebody might come up with a terrific idea.
Thanks for everyone being involved so far!
Cheers,
PI
-
- Posts: 105
- Joined: Mon Jul 27, 2009 4:06 pm
- Location: Cambridge, MA
This is a constrained optimization problem. I'm quite impressed with how it didn't work out to a closed form solution, though disappointed in the amount of time it took me to figure that out. I've derived the solution for you using two different numerical methods and provided matlab scripts. Sorry, I spent enough time on this already that I don't want to do the C++ for you.
There isn't enough space on this forum to post the results so check them out here: http://cheshirekow.com/blog/?p=86. Here is a nice picture though to get you enticed:
And a detail:
I'll do a latex and post the PDF if I get a chance tomorrow.
P.S. A very interesting problem, though the framework was easy enough, there were definitely some snags along the way (i.e. trying to factor that polynomial).
There isn't enough space on this forum to post the results so check them out here: http://cheshirekow.com/blog/?p=86. Here is a nice picture though to get you enticed:
And a detail:
I'll do a latex and post the PDF if I get a chance tomorrow.
P.S. A very interesting problem, though the framework was easy enough, there were definitely some snags along the way (i.e. trying to factor that polynomial).
Hey, that was cool :-) Nice work.
IRC: #irrlicht on irc.libera.chat
Code snippet repository: https://github.com/mzeilfelder/irr-playground-micha
Free racer made with Irrlicht: http://www.irrgheist.com/hcraftsource.htm
Code snippet repository: https://github.com/mzeilfelder/irr-playground-micha
Free racer made with Irrlicht: http://www.irrgheist.com/hcraftsource.htm
Hi! I've got some problem with the code. Please note that this is only a quick-and-dirty conversion from your matlab script into cpp.
I've tried the steepest descent method first. I've maximized the number of iterations, just in case, and it frequently hits the limit. Plus, no matter how many iterations it took to compute, it'll always return the intersection point or an other one in the extreme, as the closest point.
I've marked the problematic part with stars. I guess you're getting the imaginary part of q.Z, but I haven't treated q.Z as a complex number. Is there a workaround for it or I have to use complex numbers?
Or can you see anything else that might be an issue?
Cheers,
PI
I've tried the steepest descent method first. I've maximized the number of iterations, just in case, and it frequently hits the limit. Plus, no matter how many iterations it took to compute, it'll always return the intersection point or an other one in the extreme, as the closest point.
I've marked the problematic part with stars. I guess you're getting the imaginary part of q.Z, but I haven't treated q.Z as a complex number. Is there a workaround for it or I have to use complex numbers?
Or can you see anything else that might be an issue?
Code: Select all
core::vector3df CGame::getClosestPointToEllipsoid(const core::vector3df& Point, const core::vector3df& Scale, int &iterations)
{
// create an initial guess by finding the intersection of the ray from the
// centroid to point with the ellipsoid surface
float a = Scale.X;
float b = Scale.Y;
float c = Scale.Z;
core::vector3df p = Point;
float k = sqrt( 1 / ( (p.X*p.X)/(a*a) + (p.Y*p.Y)/(b*b) + (p.Z*p.Z)/(c*c) ) );
core::vector3df q = k * p;
// this is the scale factor for our movement along the ellipse, this can
// start off quite large because it will be reduced as needed
float s = 1.f;
// this is the main loop, you'll want to put some kind of tolerance based
// terminal condition, i.e. when the cost function (distance) is only
// bouncing back and forth between a tolerance or something more tailored to
// your application, I will stop when the distance is +/- 1%
core::vector3df plotQ[100];
float plotJ[100];
bool quit = false;
int i = 1;
do
{
// calculate the z for our given x and y
float qz_plus = c * (float)sqrt( 1.f - ((q.X*q.X)/(a*a)) - ((q.Y*q.Y)/(b*b)) );
float qz_minus = -qz_plus;
// we want the one that get's us closest to the target so
if ( fabs(p.Z - qz_plus) < fabs(p.Z - qz_minus) )
q.Z = qz_plus;
else
q.Z = qz_minus;
// calculate the current value of the cost function
float J = (float)sqrt( ((p.X - q.X)*(p.X - q.X)) + ((p.Y - q.Y)*(p.Y - q.Y)) + ((p.Z - q.Z)*(p.Z - q.Z)) );
// store the current values for the plots
plotQ[i] = q;
plotJ[i] = J;
// check to see if we overshot the goal or jumped off the ellipsoid
if ( i > 1 )
{
// if we jumped off the ellipsoid or overshot the minimal cost
//if( imag(qz) ~= 0 || J > Jplot(i-1) ) ****************************
if ( J > plotJ[i-1] )
{
// then go back to the previous position and use a finer
// step size
q = plotQ[i-1];
J = plotJ[i-1];
s /= 10.f;
--i;
// if we did just jump over the actual minimum, let's check to
// see how confident we are at this point, if we're with in
// 1% there's no need to continue
if ( i > 3 && fabs( (plotJ[i-1] - plotJ[i-3]) / plotJ[i-3] ) < 0.001f ) quit = true;
}
}
if (!quit)
{
// calculate the gradient of the cost with respect to our control
// variables; since we divide by qz in the second term we first need
// to check whether or not qz is too close to zero; if it is, we know
// that the second term should be zero
float dJdx, dJdy;
//if( q.Z < 1e-10 )
if (core::equals(q.Z, 0.f))
{
dJdx = -2.f*(p.X - q.X);
dJdy = -2.f*(p.Y - q.Y);
}
else
{
dJdx = -2.f*(p.X - q.X) + 2*(q.Z < 0.f ? -1.f : 1.f)*(p.Z - q.Z)*( c/(a*a) * (q.X/q.Z) );
dJdy = -2.f*(p.Y - q.Y) + 2*(q.Z < 0.f ? -1.f : 1.f)*(p.Z - q.Z)*( c/(b*b) * (q.Y/q.Z) );
}
// calculate the update vector that we will move along
float dx = -J/dJdx, dy = -J/dJdy;
// calculate the magnitude of that update vector
float magnitude = (float)sqrt( dx*dx + dy*dy );
// normalize our update vector so that we don't shoot off at an
// uncontrolable rate
dx = s * (1.f/magnitude) *dx;
dy = s * (1.f/magnitude) *dy;
// update the current position
q.X = q.X + dx;
q.Y = q.Y + dy;
}
// increment the index
if (++i >= 100) quit = true;
}
while (!quit);
iterations = i;
return q;
}
PI
-
- Posts: 105
- Joined: Mon Jul 27, 2009 4:06 pm
- Location: Cambridge, MA
From first glance, I don't see anything else that could be an issue, it seems you recreated the code faithfully.
You don't have to use complex numbers, as the only important thing is that you have to handle the special case if z is complex, but you don't actually have to use the value. Just check to see if the argument of the sqrt() is negative, and if it is you know you'll have a complex z and you've jumped off the ellipse.
Also, if you're point is very close to x-y plane (very low z value) when the ellipse is centered at the origin, using z as the dependent variable might not be a very good choice because you'll end up jumping off the ellipse a lot. I actually edited my article to add that note. For your test, make sure that you don't have a small z-value in your target point. This shouldn't prevent finding the solution, but it will lead to a larger number of iterations, as it will take a much finer degree of movement in the x,y directions to get that last little bit of the cost down.
You don't have to use complex numbers, as the only important thing is that you have to handle the special case if z is complex, but you don't actually have to use the value. Just check to see if the argument of the sqrt() is negative, and if it is you know you'll have a complex z and you've jumped off the ellipse.
Also, if you're point is very close to x-y plane (very low z value) when the ellipse is centered at the origin, using z as the dependent variable might not be a very good choice because you'll end up jumping off the ellipse a lot. I actually edited my article to add that note. For your test, make sure that you don't have a small z-value in your target point. This shouldn't prevent finding the solution, but it will lead to a larger number of iterations, as it will take a much finer degree of movement in the x,y directions to get that last little bit of the cost down.