Raytracing problem - Nearest point on ray
Author
I''m a somewhat experienced programmer, trying to test my wits by writing up a small and dirty raytracer with as few external references as possible, and here''s my first obstacle...
If I have a point (call it the center of an object) and a ray described by a start point and a direction vector (or by two points on an infinite line, I''ll just ignore behind-origin hits), how do I find the point on that ray/line which is closest to the center point I''ve got?
There''s lots of tutorials around describing how to find the nearest distance, but none I can find describing how to find the position of the nearest point itself.
I''d prefer practical explanations over straight algebra, please.
>>I'd prefer practical explanations over straight algebra, please.
Sorry, I barely understand how this works (or if it works)..
[edited by - mattbbangin on February 28, 2004 7:00:04 AM]
Sorry, I barely understand how this works (or if it works)..
static inline Vector3f ISOMath_RayfClosestPointToVector3f(Rayf l, Vector3f p){ Vector3f c; c.x = p.x-l.start.x; c.y = p.y-l.start.y; c.z = p.z-l.start.z; Vector3f V; V.x = l.finish.x-l.start.x; V.y = l.finish.y-l.start.y; V.z = l.finish.z-l.start.z; float d = ISOMath_Vector3fLength(V); V = ISOMath_Vector3fNormalize(V); float t = ISOMath_Vector3fDotProduct(V,c); if(t<0.0) return l.start; if(t>d) return l.finish; V.x = V.x*t; V.y = V.y*t; V.z = V.z*t; Vector3f result; ADD_VEC(l.start,V,result); return(result);}[edited by - mattbbangin on February 28, 2004 7:00:04 AM]
Author
That code certainly doesn't work: pumping in the ray {(-10,5,0) to (10,5,0)} and the point (0,0,0) returns (9,5,0), when the obvious real answer is (0,5,0).
Unless it's my dot product code that doesn't work...
[edited by - Wyrframe on February 28, 2004 2:37:25 PM]
Unless it's my dot product code that doesn't work...
float Vector::operator % (Vector &other) const { float rv; rv = (x * other.x) + (y * other.y) + (z * other.z); rv = (float)( acos(rv) ); return rv;}[edited by - Wyrframe on February 28, 2004 2:37:25 PM]
I don''t really know the actual algebra involved with this, but I think I might be able to describe sorta how to do it...
So you have a line segment (L) in space, and a point (P) somewhere else... The slope of L is m. The slope of a line that is normal to the line is m'' (-1/m). If you form a line that passes through P with slope m'', then you just need to find an intersection between the two lines.
The only thing then is the ends of the line segment. If the intersection occurs outside of the line segment, then just clamp the intersection to the endpoint.
That should work... haven''t actually gone and tested it, but it looks about right.
So you have a line segment (L) in space, and a point (P) somewhere else... The slope of L is m. The slope of a line that is normal to the line is m'' (-1/m). If you form a line that passes through P with slope m'', then you just need to find an intersection between the two lines.
The only thing then is the ends of the line segment. If the intersection occurs outside of the line segment, then just clamp the intersection to the endpoint.
That should work... haven''t actually gone and tested it, but it looks about right.
p1: vector, initial point of ray
d: vector, unit direction of ray
t: scalar, on [0, infinity) just a parameter for the equation
x: vector, the point of interest
s: scalar, distance from any point on the ray to x
p(t) = p1 + t*d // p(t) is any point on the ray
s(t) = |x - p(t)| // s(t) is the distance from x to p(t)
From here all you want to do is find t such that s(t) is at its smallest value. If you graph s as a function of t you will find that it has one point along the t axis where ds/dt is zero. This is your min point (s is minimized). You can find this point by...
Taking the derivative of s with respect to t, setting it equal to zero and solving for t (call the answer t'').
The point on the ray that is closest to x is then p(t'').
d: vector, unit direction of ray
t: scalar, on [0, infinity) just a parameter for the equation
x: vector, the point of interest
s: scalar, distance from any point on the ray to x
p(t) = p1 + t*d // p(t) is any point on the ray
s(t) = |x - p(t)| // s(t) is the distance from x to p(t)
From here all you want to do is find t such that s(t) is at its smallest value. If you graph s as a function of t you will find that it has one point along the t axis where ds/dt is zero. This is your min point (s is minimized). You can find this point by...
Taking the derivative of s with respect to t, setting it equal to zero and solving for t (call the answer t'').
The point on the ray that is closest to x is then p(t'').
hmmm, that derivative is non-trivial.
What BlabberBoy is talking about will sort of work. You have to keep in mind that lines in R3 don''t really have slopes like lines in R2. So finding a line that passes through your ''center point'' and is perpendicular to your ray is going to be tricky.
One way to do this...
k1 = d cross (p1-x)
normalize(k1)
k2 = k1 cross d
k2 is a unit vector that is perpendicular to your ray and points toward x.
Now you have two ray equations:
p(t) = p1 + t*d
o(t) = x + t*k1
Set them equal to each other, then solve for t, again call the answer t''. p(t'') is the point along your ray that is closest to x.
What BlabberBoy is talking about will sort of work. You have to keep in mind that lines in R3 don''t really have slopes like lines in R2. So finding a line that passes through your ''center point'' and is perpendicular to your ray is going to be tricky.
One way to do this...
k1 = d cross (p1-x)
normalize(k1)
k2 = k1 cross d
k2 is a unit vector that is perpendicular to your ray and points toward x.
Now you have two ray equations:
p(t) = p1 + t*d
o(t) = x + t*k1
Set them equal to each other, then solve for t, again call the answer t''. p(t'') is the point along your ray that is closest to x.
Author
All would be fine and good if it weren''t for the fact that computers cannot "set equal to zero and solve". And nobody has yet taught me how to re-write that into "add this, multiply this, and there you go".
Computers can do nothing but straight arithmetic, and so that''s the realm that one must work in when working with them. And besides which, you''re looking at an absolute value; the point we''re looking for on the derivative of s() has no unique solution at the sharp point.
And I know for certain now that my dot product is a bit messed. Could somebody point me to a decent implementation? I can''t find any definition that doesn''t require you to already know the angle between them.
Computers can do nothing but straight arithmetic, and so that''s the realm that one must work in when working with them. And besides which, you''re looking at an absolute value; the point we''re looking for on the derivative of s() has no unique solution at the sharp point.
And I know for certain now that my dot product is a bit messed. Could somebody point me to a decent implementation? I can''t find any definition that doesn''t require you to already know the angle between them.
a dot b = a.x*b.x + a.y*b.y + a.z*b.z
[edited by - thedustbustr on February 28, 2004 4:55:34 PM]
float DotProduct(const vector3 &a, const vector3 &b) { return a.x*b.x + a.y*b.y + a.z*b.z; } [edited by - thedustbustr on February 28, 2004 4:55:34 PM]
This problem is simple enough. If the ray starts at point O with direction d and we want to find the point of the ray that is closest to P, you just have to compute O+((P-O).d)*d/Norm(d). In C++:
[edited by - Alvaro on February 28, 2004 4:35:16 PM]
#include <iostream>struct Vector { double x,y,z; Vector(double x,double y,double z):x(x),y(y),z(z){}};Vector operator+(Vector const &v1, Vector const &v2){ return Vector(v1.x+v2.x,v1.y+v2.y,v1.z+v2.z);}Vector operator*(double s, Vector const &v){ return Vector(s*v.x,s*v.y,s*v.z);}Vector operator*(Vector const &v, double s){ return s*v;}Vector operator/(Vector const &v, double s){ double i=1.0/s; return i*v;}double operator*(Vector const &v1, Vector const &v2){ return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;}double Norm(Vector const &v){ return v*v;}struct Point { double x,y,z; Point(double x,double y,double z):x(x),y(y),z(z){}};std::ostream &operator<<(std::ostream &o, Point const &p){ return o << '(' << p.x << ',' << p.y << ',' << p.z << ')';}Vector operator-(Point const &p1, Point const &p2){ return Vector(p1.x-p2.x,p1.y-p2.y,p1.z-p2.z);}Point operator+(Point const &p, Vector const &v){ return Point(p.x+v.x,p.y+v.y,p.z+v.z);}struct Ray { Point p; Vector v; Ray(Point const &p, Vector const &b):p(p),v(v){} Ray(Point const &p1, Point const &p2):p(p1),v(p2-p1){}};Point closest(Ray const &r, Point const &p){ return r.p+(p-r.p)*r.v*r.v/Norm(r.v);}int main(void){ std::cout << closest(Ray(Point(-10,5,0),Point(10,5,0)),Point(0,0,0)) << std::endl; return 0;}[edited by - Alvaro on February 28, 2004 4:35:16 PM]
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
