Author Topic: Ray intersect Sphere?? (Read 10182 times)

Devils Child · « **on:** August 07, 2007 »

hi there...
i'm thinking of making a raytracer again. this time, a better one. so i need a triangle-sphere intersection...

where is one?
here it is: http://www.blitzbasic.com/codearcs/codearcs.php?code=947

well then, i implemented this one pretty well, but i have 2 problems:

1. ==== backward intersection ====
the line which intersects the sphere, hits the sphere even when the sphere is "behind" the line. it looks like "backward projection" then...
2. ==== the reflectional ray ====
what i need to calculate is a ray that is reflected by the sphere. this reflectional ray starts at the intersection point between the line and the sphere and the direction it goes, is what i need to calculate.

if you would be so good and help me

i would really appreciate any help coming here

edit: my try
(original code by elias_t)

Code: [Select]

;infinite ray to sphere intersection example.
;The intersection points are returned.

Dim picked1#(2);1st intersection point
Dim picked2#(2) ;<=== the reflection vector!!!


; x1,y1,z1    : 1st point of segment
; x2,y2,z2    : 2nd point of segment
; x3,y3,z3, r : coordinates And radius of sphere)

Function sphere_line_intersection (x1#,y1#,z1#, x2#,y2#,z2#, x3#,y3#,z3#,r#)

Local mu#

Local dx1#=(x2-x1)
Local dy1#=(y2-y1)
Local dz1#=(z2-z1)

Local dx2#=(x1-x3)
Local dy2#=(y1-y3)
Local dz2#=(z1-z3)


Local a# = dx1*dx1 + dy1*dy1 + dz1*dz1
Local b# = 2*( dx1*dx2 + dy1*dy2 + dz1*dz2 )
Local c# = (x3*x3)+(y3*y3)+(z3*z3)+(x1*x1)+(y1*y1)+(z1*z1)-2*(x3*x1+y3*y1+z3*z1)-(r*r)


Local bm#=(b * b)
Local da#=2*a
Local fac#=(4 * a * c)

Local i#=bm - fac

Local sqi#=Sqr(i)



;if segment crosses the sphere
If i>=0

  ;incoming intersection
  mu = (-b - sqi) / da
  picked1(0) = x1 + mu*dx1
  picked1(1) = y1 + mu*dy1
  picked1(2) = z1 + mu*dz1



  ;============================
  ;====== reflection ray ======
  ;============================
  <=== the reflection vector!!!
  
  mu = (-b + sqi) / da
  picked2(0) =  x1# + picked1(0)
  picked2(1) =  y1# + picked1(1)
  picked2(2) =  z1# + picked1(2)
  ;

Return 1
EndIf



End Function







;************************************************************************
;EXAMPLE




Graphics3D 640,480,0,2

cam=CreateCamera()
MoveEntity cam,0,3,-5
light=CreateLight()


;------------------------------
;create a "ray like triangle" near the camera's position 0,0,-5
ray=CreateMesh()
rs=CreateSurface(ray)
v0=AddVertex(rs,-.03,0,-5)
v1=AddVertex(rs,0,0,10)
v2=AddVertex(rs,.03,0,-5)
AddTriangle(rs,v0,v1,v2)
EntityColor ray,255,255,0
UpdateNormals ray
EntityFX ray,17
;------------------------------


;------------------------------
ray2=CreateMesh()
rs2=CreateSurface(ray2)
v02=AddVertex(rs2,-.03,0,-5)
v12=AddVertex(rs2,0,0,+10)
v22=AddVertex(rs2,.03,0,-5)
AddTriangle(rs2,v02,v12,v22)
EntityColor ray2,0,255,0
UpdateNormals ray2
EntityFX ray2,17
;------------------------------


;------------------------------
;Create a sphere [radius=1]
ob1=CreateSphere()
EntityColor ob1,255,0,0
EntityAlpha ob1,.25
EntityFX ob1,16
ScaleEntity ob1,1.5,1.5,1.5

;create dummy objects
d1=CreateCube()
ScaleEntity d1,.05,.05,.05
EntityColor d1,0,0,255
d2=CreateCube()
ScaleEntity d2,.05,.05,.05
EntityColor d2,0,255,0
;------------------------------


;variables to control the ray
rx#=0
ry#=0
rz#=10


PointEntity cam,ob1

	
While Not KeyHit(1)
RenderWorld()

;move the ray
If KeyDown(200) ry=ry+.1 VertexCoords rs,v1,rx,ry,rz
If KeyDown(208) ry=ry-.1 VertexCoords rs,v1,rx,ry,rz
If KeyDown(203) rx=rx-.1 VertexCoords rs,v1,rx,ry,rz
If KeyDown(205) rx=rx+.1 VertexCoords rs,v1,rx,ry,rz

VertexCoords rs2, v02, picked1(0) - .03, picked1(1), picked1(2)  ;=AddVertex(rs2,-.03,0,-5)
VertexCoords rs2, v12, picked1(0) + picked2(0) - .03, picked1(1) + picked2(1), picked1(2) + picked2(2)
VertexCoords rs2, v22,  picked1(0) + .03, picked1(1), picked1(2) ;=AddVertex(rs2,.03,0,-5)

;x,y,z=1st point ,rx,ry,rz=2nd point ,sx,sy,sz,r=sphere pos and radius r
If sphere_line_intersection (0,0,-5, rx,ry,rz, 0,0,0,1.5)
PositionEntity d1,picked1(0),picked1(1),picked1(2)
PositionEntity d2,picked2(0),picked2(1),picked2(2)
EndIf

Text 200,0,"[ Cursor keys change ray direction]"
Text 200, 20, "The yellow ray is the ray."
Text 200, 40, "The green ray is the reflective ray."
Flip
Wend
ClearWorld()
End

Stonemonkey · « **Reply #1 on:** August 07, 2007 »

Sorry, I'm not sure I understand the maths in that but I'm pretty sure the reflected vector shouldn't be calculated using x1,y1,z1.

Here's an attempt using a different method (probably a bit slower due to quite a few SQRs)but I'm not sure it's working as you want as the green cube doesn't seem to be in the right place.

Code: [Select]

Function sphere_line_intersection2 (x0#,y0#,z0#, x1#,y1#,z1#,cx#,cy#,cz#,r#)

;get ray vector and normalise
	nx#=x1-x0
	ny#=y1-y0
	nz#=z1-z0
	d#=1.0/Sqr(nx*nx+ny*ny+nz*nz)
	nx=nx*d
	ny=ny*d
	nz=nz*d
	
;get vector from one point on the ray to the centre of the sphere
	vx#=cx-x0
	vy#=cy-y0
	vz#=cz-z0
	
;dot product to find distance from point 0 to the point of the line closest to the centre of the sphere
	d#=vx*nx+vy*ny+vz*nz
	
;is the point in front?
	If d>0.0 Then
	
;find the coordinates of the closest point of the ray to the centre of the sphere
		px#=x0+nx*d
		py#=y0+ny*d
		pz#=z0+nz*d
	
;find the distance from the centre of the sphere to the closest point on the ray
		d#=Sqr((px-cx)^2+(py-cy)^2+(pz-cz)^2)
	
;does the ray intersect the sphere?
		If d<=r Then
		
;use pythagoras to find 1/2 * the size of the chord (how much of the ray is within the sphere)
			d#=Sqr(r*r-d*d)
		
;subtract that from the closest point (back towards point 0)
			picked1(0)=px-nx*d
			picked1(1)=py-ny*d
			picked1(2)=pz-nz*d
		
;find the vector from the centre of the sphere to the point of intersection and normalise (the collision normal) 
			vx=picked1(0)-cx
			vy=picked1(1)-cy
			vz=picked1(2)-cz
			d=1.0/Sqr(vx*vx+vy*vy+vz*vz)
			vx=vx*d
			vy=vy*d
			vz=vz*d
			
;reflect the ray vector n using collision normal v
			d=2.0*(vx*nx+vy*ny+vz*nz)
			picked2(0)=nx-vx*d
			picked2(1)=ny-vy*d
			picked2(2)=nz-vz*d
			
			Return 1
			
		End If
		
	End If
	
	Return 0
	
End Function

Stonemonkey · « **Reply #2 on:** August 07, 2007 »

About the problem with the green cube, maybe multiplying the reflected vector by some amount makes it look like it should.

replace the 3 lines in the function above with:

Code: [Select]

			picked2(0)=(nx-vx*d)*50.0
			picked2(1)=(ny-vy*d)*50.0
			picked2(2)=(nz-vz*d)*50.0

Devils Child · « **Reply #3 on:** August 07, 2007 »

holy jesus that is cool!
it works. it is 30% slower, but works

http://www.dc.freecoder-portal.de/upload/files/Devils%20Child/Screen190.jpg

and it is not picking backward.

the only thing left to do is:
- infinite ray
- and refraction

thank you so much for this code

i love you!

Stonemonkey · « **Reply #4 on:** August 07, 2007 »

If you're planning on several reflections per ray then you'd only need to normalise the initial ray before calling the function the first time instead of normalising at the start of the function each time as the reflected ray is normalised anyway.

There should be some more optimisations that can be done too but I left it at that as it's a bit clearer on how it works.

Stonemonkey · « **Reply #5 on:** August 07, 2007 »

This might be a little better, the normalisation at the start of the function can still be removed and only done for the initial raycast but I've managed to remove 2 of the other SQRs so it only needs 1 per reflected ray And you could even get rid of the divide by r and multiply by 2.0 in the reflection calcs by storing 2.0/r along with the sphere radius.

Code: [Select]

Function sphere_line_intersection2 (x0#,y0#,z0#, x1#,y1#,z1#,cx#,cy#,cz#,r#)

;get ray vector and normalise
	nx#=x1-x0
	ny#=y1-y0
	nz#=z1-z0
	d#=1.0/Sqr(nx*nx+ny*ny+nz*nz)
	nx=nx*d
	ny=ny*d
	nz=nz*d

;get vector from one point on the ray to the centre of the sphere
	vx#=cx-x0
	vy#=cy-y0
	vz#=cz-z0

;dot product to find distance from point 0 to the point of the line closest to the centre of the sphere
	d#=vx*nx+vy*ny+vz*nz

;is the point in front?
	If d>0.0 Then

;find the coordinates of the closest point of the ray to the centre of the sphere
		px#=x0+nx*d
		py#=y0+ny*d
		pz#=z0+nz*d

;find the square of the distance from the centre of the sphere to the closest point on the ray
		d#=(px-cx)*(px-cx)+(py-cy)*(py-cy)+(pz-cz)*(pz-cz)

;does the ray intersect the sphere?
		If d<=r*r Then
		
;use pythagoras to find 1/2 * the size of the chord (how much of the ray is within the sphere)
			d#=Sqr(r*r-d)

;subtract that from the closest point (back towards point 0)
			picked1(0)=px-nx*d
			picked1(1)=py-ny*d
			picked1(2)=pz-nz*d

;find the vector from the centre of the sphere to the point of intersection no need to normalise as it's length is equal to r (may need to normalise for lighting) 
			vx=picked1(0)-cx
			vy=picked1(1)-cy
			vz=picked1(2)-cz

;reflect the ray vector n using collision normal v
			d=2.0*(vx*nx+vy*ny+vz*nz)/r
			picked2(0)=nx-vx*d
			picked2(1)=ny-vy*d
			picked2(2)=nz-vz*d

			Return 1

		End If

	End If

	Return 0

End Function

Devils Child · « **Reply #6 on:** August 07, 2007 »

while i have your attention: can i ask you some questions please?

my first question is how to get lighting work.
in my raytracer i have several point lights hanging arround somewhere, so how can i calculate the intensity of a light to a normal/position of the sphere?

i mean - i only have the position, the normal of this position and the light position given. i need to calculate the light intensity there (including shininess if possible)

thanks a lot

edit: indeed, your update is down from 1500 millisecs to 1000 millisecs rendertime! faster!

Stonemonkey · « **Reply #7 on:** August 07, 2007 »

The way I usually do it, I'm not saying this is all correct and there are probably many ways to do it and I might have some of the terms wrong.

Diffuse lighting: I use (normalised vector from point->light)dot(surface normal) which will return a value -1.0 to 1.0

I either clip off the negative part (If diffuse<0.0 then diffuse=0.0)but that can result in very dark/black patches or scale it and move it into the range 0.0 to 1.0 using(diffuse=0.5+diffuse*0.5)
and the final colour may be sphere_colour*light_colour*diffuse (I have colour components in the range 0.0 to 1.0 when I do this)

Specular lighting: (normalised vector from point->light)dot(reflected vector) again a value -1.0 to 1.0

Then clip it if it's negative else multiply it by itself several times, the more times the shinier the object.
and the final colour may just use the light colour (could use the object colour depending on the properties you want) light_colour*specular which is added to the diffuse final colour.

You also need a control value that starts at 1.0 and decreases somehow with each reflection relative to how reflective the object is, the more reflective the object, the less it decreases by. So again you could give a sphere a reflective value from 0.0 to 1.0 (0.0=non reflective 1.0=totally reflective) and you multiply the control value by that for each reflection and you multiply your lit colour values by the control value before adding them to the total colour for that pixel.

Sorry i'm rambling on a bit, any questions just ask.

EDIT:
And you have to test for shadows which should just be a quick intersection test of all the objects between the point and the lightsource which has just made me re think about the scaling for diffuse lighting, might not be a good idea after all as shadows from objects could appear on the wrong side of objects.

taj · « **Reply #8 on:** August 08, 2007 »

Stonemonkey,

the dark effect is probably because you dont seem to be adding in any ambient light factor (0.1 say). Rather than scale the diffuse to fit I'd usually go for:

Ka + Kd*clamp(diffuse,0,1) ...plus specular

Also if you do this per light or per scene it makes a difference. If you add in Ka per light, scenes can get washed out by a lot of lights but maybe thats what we would want? Theoretically, as Ka is a basic approximation to the diffuse interreflections in the scene, each new light does add in more.

---- theory ramblings

So raytracing doesnt calculate diffuse interreflections like say radiosity. But I had an idea this morning to overcome this partially. The back side of a surface is going to get no direct light (diffuse and specular) only ambient. The question in my mind was can ambient be calculated more accurately than just Ka? Well given its the diffuse interreflection where is my sphere going to get most indirect light from - which direction? From another surface directly facing the light (hence most bright). Where are those surfaces relative to my sphere? They are in the opposite direction of the light from my sphere normal. ie the most ambient would be where the diffuse component of my sphere is -1.0!!! My sphere is facing away from the light so its "looking at" surfaces facing the light! Similarly, Least ambient would be where diffuse = 1.0 (ie the diffuse surfaces it faces at that point would be facing away from the light).

Sooo... a better approximation for ambient lighting would be:

Ka * (1-(diffuse+1)*0.5) + diffuse + specular (for each light)

This gives no ambient lighting contribution in the specular area of the sphere and maximum at the opposite side but it generalises to any surface of course. This is actually quite close to what your are doing with your diffuse lighting Stonemonnkey but I think its more mathematically correct. Would need some fine tuning I guess.

Did I manage to explain it? Its sort of the opposite of ambient occlusion I think and could be used with it. No colour bleeding :-( like radiosity.

Chris

Devils Child · « **Reply #9 on:** August 08, 2007 »

for the lighing question:

i know how shadows work, but first i need lighting.

this screenshot should explain what i ask for:
http://dc.freecoder-portal.de/upload/files/Devils%20Child/Screen191.jpg

Stonemonkey · « **Reply #10 on:** August 08, 2007 »

For diffuse, calculate the dot product of that vector (unit vector pointing towards the light) and the surface normal (unit vector pointing directly away from the centre of the circle)

For specular, calculate the dot product of that vector (unit vector pointing towards the light) and the reflected vector.

Devils Child · « **Reply #11 on:** August 08, 2007 »

excuse me: how do i do a dot product again?

Stonemonkey · « **Reply #12 on:** August 08, 2007 »

There's reference to it in the function but here you go anyway:

Code: [Select]

dot_product=vx*nx+vy*ny+vz*nz

if both vectors are unit vectors then the result is in the range -1.0 to 1.0 and is the cosine of the angle between them

so:

Code: [Select]

angle=acos(vx*nx+vy*ny+vz*nz)

but you'll rarely need to find the angle.

Devils Child · « **Reply #13 on:** August 08, 2007 »

ok. got that understand.
i will implement it right now!

this is a screen so far:
http://www.dc.freecoder-portal.de/upload/files/Devils%20Child/Screen196.jpg

Jim · « **Reply #14 on:** August 09, 2007 »

K+ for Stonemonkey - :goodpost:s!

Devils Child · « **Reply #15 on:** August 09, 2007 »

i agree with that: karma for Stonemonkey!

am i pushy when i ask how i can transform this function to a cube or a cylinder?

Stonemonkey · « **Reply #16 on:** August 09, 2007 »

Yes. Just kidding, no but it's not really something I've looked into myself much. I did have a go at cylinders once but never quite worked it out. You've given me the raytracing bug a bit now and I'm having a go myself so I'll have a think about that and see if i can come up with anything.

Devils Child · « **Reply #17 on:** August 10, 2007 »

cool! so you like raytracing, hm?
what i still need is a cube and a cylinder.
i suppose a tube or a torus would be hard

Stonemonkey · « **Reply #18 on:** August 10, 2007 »

for a tube or cylinder you'll use the vector running up the centre of it and find if the ray passes that line closer than the radius so it'll involve finding the closest points of 2 vectors, if they're closer than radius then you have to find the point on the ray that the distance to the closest point on the cylinder axis is equal to radius then find if that is between the ends of the cylinder.

for a torus you have a circle and you have to do a similar thing with the edge of that circle

putting it all into practice is a different matter though.

Devils Child · « **Reply #19 on:** August 10, 2007 »

but a cube could be pretty easy tho?
btw: i still dont understand picking algorithms in general. i'm in class 10!! we didnt made picking algorithms in school, but still, i wonder why i dont understand these...

i found a cube picking algorith here:
http://www.blitzbasic.com/codearcs/codearcs.php?code=1029

the first matter is that there are no picked coordinate returns. also a reflection ray would be useful. yet, i dont need reflection. i have other problems at the moment.

basicaly what i need to do is to make this cube picking algorithm go:
- with infinite picking line
- NO backward picking (you know what i mean...)
- reflection (i will need that later)

btw: i really appreciate that you help me a lot!
i like this community somehow... k+ for stonemonkey!