Sign in to follow this  
ToohrVyk

[Math] Ray-Sphere

Recommended Posts

I wish to determine the parameters tin and tout where a ray enters and leaves a sphere. Is there a formula to obtain these values without computing a square root?

Share this post


Link to post
Share on other sites
no. If there would be, then we would be able to use it to compute square root without computing square root. [lol]

btw if you need code,

//
// Line-sphere intersection.
// p=(line origin position - sphere position),
// d=line direction,
// r=sphere radius,
// i1=first intersection distance,
// i2=second intersection distance
// i1<=i2
//
// returns true if intersection found,false otherwise.
//
bool LineSphereIntersect(const Vec3d &p, const Vec3d &d,real64 r, real64 &i1, real64 &i2){
real64 det,b;
b = -DotProd(p,d);
det=(b*b) - (p.x*p.x + p.y*p.y + p.z*p.z) + r*r;
if (det<0){
return false;
}//else
det= sqrt(det);
i1= b - det;
i2= b + det;
return true;

};





(and of course you can compute that not using sqrt and using something else instead, like det=exp(log(det)/2.0)) but it's square root anyway)
As you can see from formula , the "det" is kind of squared distance from middle between intersection points, so maybe you can still avoid square root somehow.

edit: By the way, it doesn't compute square root for rays that do not intersect sphere, so if most rays do not intersect your spheres it is almost as good as if it wouldn't compute square root at all.

Share this post


Link to post
Share on other sites
Quote:
Original post by Dmytry
no. If there would be, then we would be able to use it to compute square root without computing square root.[lol]


Excellent answer [lol]

Quote:
btw if you need code
*** Source Snippet Removed ***


Thank you, but your code does not apply for three reasons: I'm not using a C family language, my ray direction vector is not normalized, and I'm actually solving an equivalent problem (moving sphere vs moving sphere collision).

Snippet of my (untested) code:

(** Compute the two tangential contact times for a two-sphere
system.
@param pA The initial (x,y) position of the first sphere.
@param vA The velocity of the first sphere.
@param rA The radius of the first sphere.
@param pB The position of the second sphere.
@param vB The velocity of the second sphere.
@param rB The radius of the second sphere.
@return The times (t0,t1) of tangential contact.
@raise NoCollision If there is no contact. *)
let sphereContact (pA) (vA) (rA) (pB) (vB) (rB) =
let px = (fst pA) -. (fst pB)
and py = (snd pA) -. (snd pB) in
let vx = (fst vA) -. (fst vB)
and vy = (snd vA) -. (snd vB) in
let r = rA +. rB in
let sq x = x *. x in

let delta =
let sub = px *. vy -. py *. vx in
(sq r) *. ((sq px) +. (sq py)) -. (sq sub)
in

if delta < 0. then raise NoCollision;

let d = ((sq vx) +. (sq vy)) in
let delta = (sqrt delta) /. d in
let x = -. (px *. vx +. py *. vy) /. d in

(x -. delta),(x +. delta);;

Share this post


Link to post
Share on other sites
Ray-sphere collision detection is pretty easy and can be found here. Moving sphere-sphere collision isn't that bad either and can be found here. The latter link only checks to see if two spheres will intersect, but it provides you with the second-order equation you would solve to get intersection times so you can modify the algorithm.

Share this post


Link to post
Share on other sites
Quote:
Original post by ToohrVyk
Quote:
Original post by Dmytry
no. If there would be, then we would be able to use it to compute square root without computing square root.[lol]


Excellent answer [lol]

Quote:
btw if you need code
*** Source Snippet Removed ***


Thank you, but your code does not apply for three reasons: I'm not using a C family language, my ray direction vector is not normalized,

modded for non-normalized direction vector:


//
// Line-sphere intersection.
// p=(line origin position - sphere position),
// d=line direction,
// (line equation is p+d*t)
// r=sphere radius,
// i1=first intersection,
// i2=second intersection
// i1<=i2
//
// returns true if intersection found,false otherwise.
//
bool LineSphereIntersect(const Vec3d &p, const Vec3d &d, real64 r, real64 &i1, real64 &i2){
real64 det,b;
b = -DotProd(p,d);
real64 scale = 1.0/DotProd(d,d);
b*=scale;
det=(b*b) - ((p.x*p.x + p.y*p.y + p.z*p.z) + r*r)*scale;
if (det<0){
return false;
}//else
det= sqrt(det);
i1= b - det;
i2= b + det;
return true;
};






Not tested, but I think should work.

Quote:


and I'm actually solving an equivalent problem (moving sphere vs moving sphere collision).

Yes, that does reduce to ray-sphere intersection for sphere with radius equal to sum of radiuses.

Quote:


Snippet of my (untested) code:

*** Source Snippet Removed ***

OCaml? something else ml-based?
no offense, but it's approximately as readable as assembler (of some procesor unknown to me). Well, i can kind of figure what is going on, but, why "+." ? Is it to show that it's floating point addition?
edit: fixed typo in code

[Edited by - Dmytry on November 27, 2005 3:48:29 PM]

Share this post


Link to post
Share on other sites
Quote:
Original post by Dmytry
OCaml? something else ml-based?
no offense, but it's approximately as readable as assembler (of some procesor unknown to me). Well, i can kind of figure what is going on, but, why "+." ? Is it to show that it's floating point addition?
edit: fixed typo in code


You're right on both accounts. It's OCaml, and +. is the floating-point variant for +

Share this post


Link to post
Share on other sites
[sorry for offtopic] Actually I just seen some sample code with Ocaml a while ago that had these strange "let" "in" "." things, but had to guess what dots everywhere is for. I have few questions: what happens if . is missed? "omg_wtf_cast" to integer? (i suppose idea is that + can't or should not be overloaded to work with integers, reals, vectors, etc.) I googled for Ocamn now, must admit some language features does look nice, but syntax looks bit like "parallel evolution" to c++ , in sense of having strange strings of symbols that is hard to decipher. (and i don't like c/c++ sintax either)

[Edited by - Dmytry on November 28, 2005 2:40:05 PM]

Share this post


Link to post
Share on other sites
OCaml is a lot more readable with decent syntax coloring and enough coding experience. And my code isn't exactly adhering to common beauty guidelines.

The difference between (+) and (+.) is their type:

val (+) : int -> int -> int
val (+.) : float -> float -> float

And since all casts in O'Caml must be explicit, the compiler will throw a fit if you use one instead of the other. There is no operator overloading in O'Caml, because the type of the operator is used to infer the type of the objects it is applied to (and not the reverse).

Here, perhaps a more readable version:

let vdot (ax,ay) (bx,by) = (ax*.bx+.ay*.by);;
let vmin (ax,ay) (bx,by) = (ax-.bx),(ay-.by);;

(** Compute the parameter of the collision between a circle
and a ray.
@param p Circle center.
@param r Circle radius.
@param a Ray origin.
@param v Ray direction (not normalized)
@return The time of first contact.
@raise NoCollision if there is no contact in the future. *)
let rayCircleContact (p) (r) (a) (v) =

let p = vmin a p in
let sq x = x *. x in

(* a = vx² + vy² *)
(* b = px.vx + py.vy --- b_h = b/2 *)
(* c = px² + py² - r² *)

let a = vdot v v in
let b_h = vdot p v in
let c = vdot p p -. (sq r) in

(* det = b² - 4ac --- det_f = det/4 *)

let det_f = (sq b_h) -. a *. c in

if det_f < 0. then raise NoCollision;

(* d = sqrt( det ) --- d_ha = d/2a *)
(* m = - b/2a *)

let d_ha = (sqrt det_f) /. a in
let m = -. b_h /. a in

(* t = - (b + d) / 2a *)

let t = (m -. d_ha) in
if t > 0. then t else raise NoCollision;;

(** Compute the first tangential contact time for a two-circle
system.
@param pA The initial (x,y) position of the first circle.
@param vA The velocity of the first circle.
@param rA The radius of the first circle.
@param pB The position of the second circle.
@param vB The velocity of the second circle.
@param rB The radius of the second circle.
@return The first time in the future when collision occurs.
@raise NoCollision If there is no contact in the future. *)
let circleCircleContact (pA) (vA) (rA) (pB) (vB) (rB) =

let p = vmin pA pB in
let v = vmin vA vB in
let r = rA +. rB in

rayCircleContact (0.0,0.0) (r) (p) (v);;

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this