Ray->Sphere intersection?

Blitz3D Forums/Blitz3D Programming/Ray->Sphere intersection?

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 realy appreciate any help comming here

edit: my try
(original code by elias_t)
;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


ummh. seems like i'm pretty much on my own right now, huh?
still waiting for someone who helps me...

What about this: http://www.blitzbasic.com/codearcs/codearcs.php?code=793

thats useless, because it intersects a polygon, not a sphere. the only matter in that code you posted is the speed

Hi Devils Child,
I've written a function to reflect a ray below. To use it you need to perform the following steps:-

1 - Create a Unit Vector called Ray from Line start x,y,z to sphere intersection point ix,iy,iz:-

    global Ray.vector=new vector
    Ray\x=ix-x : Ray\y=iy-y : Ray\z=iz-z
    UnitVector(Ray)


2 - Create a Normal Unit Vector from the sphere center sx,sy,sz to the intersection point ix,iy,iz:-

   global Normal.vector=new vector
   Normal\x=ix-sx : Normal\y=iy-sy : Normal\z=iz-sz
   UnitVector(Normal)


3 - Now you can present these to the function ReflectRay to get the new ray direction:-

  newRay.vector=ReflectRay(Normal,Ray)


Hope that helps!

Cheers,

Roy

Type vector
    Field x#,y#,z#
End Type

Function ReflectRay.vector(normal.vector,ray.vector)
	Local newray.vector=New vector
	Local normpiv%=CreatePivot()
	AlignToVector normpiv,normal\x,normal\y,normal\z,3
	TFormVector ray\x,ray\y,ray\z,0,normpiv
	ray\x=TFormedX()
	ray\y=TFormedY()
	ray\z=TFormedZ()
	TurnEntity normpiv,0,0,180
	TFormVector ray\x,ray\y,ray\z,normpiv,0
	newray\x=-TFormedX()
	newray\y=-TFormedY()
	newray\z=-TFormedZ()
	FreeEntity normpiv
	UnitVector(newray)
	Return newray
End Function

Function UnitVector(v.vector)
	Local l#
	l=Sqr(v\x*v\x+v\y*v\y+v\z*v\z)
	If l=0 Then l=1
	v\x=v\x/l
	v\y=v\y/l
	v\z=v\z/l
	If Abs(v\x)<0.000001 Then v\x=0
	If Abs(v\y)<0.000001 Then v\y=0
	If Abs(v\z)<0.000001 Then v\z=0
End Function


Solving the above for u =

(x3 - x1)(x2 - x1) + (y3 - y1)(y2 - y1) + (z3 - z1)(z2 - z1)
----------------------------------------------------------------------
(x2 - x1)(x2 - x1) + (y2 - y1)(y2 - y1) + (z2 - z1)(z2 - z1)

If u is not between 0 and 1 then the closest point is not between P1 and P2

Given u, the intersection point can be found, it must also be less than the radius r. If these two tests succeed then the earlier calculation of the actual intersection point can be applied.

robFerriby: thank you for this explaining reflection code. it explains the reflection system to me pretty well (well it's 3d code but i'll be able to do it in 2d)

wayne: thank you for this one!

besides i have asked in an other forum since i got no answer for a whole day, and they replied pretty fast and i got it working


you will hear from me when my raytracer is finished.
i know the screenie looks pretty boring. but lighting and shadowing will come!

Glad I could help. The results look pretty good so far, even at this early stage. I'll be looking forward to seeing the final version. :)

Cheers,

Roy

i think you wanna see my current one, too ;)



there are soft shadows! there are two methods to render soft shadows. the combination of both is looking very good. and it took me 18 seconds to render this.

i still need to work on the shading of the spheres, they doesnt seem to receive any light at all.

btw: the ground is NOT a plane; it's a giant sphere ;)

Do you not use clusters for raytracing? With this you will have the same quality but even with less calculations ( = faster).

See gamedev.com or any of the 3d-expert coding community for more informations. I think you increased since ~2004, so you can now realy learn from the stuff.

hm clusters? what is that?

The shadowing looks pretty cool. The spheres must be receiving some light though otherwise you'd see nothing but blackness. I noticed in your first pic that the spheres had some reflective qualities. Did you take that out?
Looking forward to the next screen shot ;)

Cheers,

Roy

i currently commented the reflection code out. i need to concentrate on the shadowing first, then i continue on reflection/refraction.


also i use a script engine for the render, so you can script your scenes such as in povRay...


edit: hm. i would really appreciate if anyone tells me where the bug is, that the spheres aren't lighted...

Remember that unit normal that you computed from the sphere center to the intersection point?
Well you'll need to perform a second unit normal which points from the intersection point to the position of the light source. Using these two normals perform a DotProduct test which should yeld a result between -1 (complete darkness) and 1 (fully lit).
Add 1 to the result and divide by 2 to give a range between 0 and 1. Now you can multiply your sphere surface RGB values against it to arrive at the correct color on the sphere surface for that particular ray position.

Cheers,

Roy

wow. i got that fixed before i read that!
will be interesting to know if i want to have shaded balls



edit: sorry, Inarie for using your scene :) its just brilliant!!

Great job! It's looking very good now :)



this is a skylight with 50x sampling

a skylight is nothing more than a light which is very big...

Looking good.

You need ray to box intersection?

This is a long shot, but I did this demo ages ago (messy), which you may be able to expand upon and use? It does ray-to-box intersection, but only ray to centre of box. Box can have any rotation and doesn't use picks.
	Graphics3D 800,600,32
	SetBuffer BackBuffer()

	WireFrame 0
	AntiAlias 0

	SeedRnd MilliSecs()

	Global frame_count%
	Global fps%
	Global slowest_fps%
	Global fps_timeout%
	Global frame_time%
	Global slowest_frame%
	Global frame_start%
	fps_timer = CreateTimer(60)
	slowmo% = False
	wiref% = False
	
	piv = CreatePivot()
	cam = CreateCamera(piv)
	PositionEntity cam,0,0,-8
	TurnEntity piv,0,45,0
	light = CreateLight()
	
	cube = CreateCube()
	EntityAlpha cube,.9
	;PositionEntity cube,-.5,-2,3
	xsize# = .5
	ysize# = 1.23
	zsize# = 1.5
	ScaleMesh cube,xsize,ysize,zsize
		
	piv2 = CreatePivot()
	enemy = CreateCone(32,1,piv2)
	RotateMesh enemy,90,0,0
	UpdateNormals enemy
	ScaleEntity enemy,.3,.3,.3
	PositionEntity enemy,1,2,-3
	EntityColor enemy,255,255,0
	
	hit = CreateSphere()
	ScaleMesh hit,.12,.12,.12
	EntityColor hit,0,255,0
	PositionEntity hit,0,0,0

	p#=-.1
	y#=-.5
	r#=1
	p2#=0
	y2#=-.5
	r2#=-1
	timer = MilliSecs()

	num_spots=100
	s=CreateSphere(3,cube)
	ScaleMesh s,.1,.1,.1
	EntityColor s,255,0,255
	HideEntity s
	Dim spots(num_spots)
	For n=1 To num_spots
		spots(n)=CopyEntity(s,cube)
		PositionEntity spots(n),0,10000,0
	Next
	free_spot=1
	
	; Main loop
	While Not KeyHit(1)
		frame_start = MilliSecs()
	
		If KeyHit(28) Then slowmo = Not slowmo
		If KeyHit(14)
			wiref = Not wiref
			WireFrame wiref
		EndIf

		PointEntity enemy,cube
		TurnEntity piv2,p2,y2,r2
		If MilliSecs()-timer>=3000
			p=Rnd(-1,1)
			y=Rnd(-1,1)
			r=Rnd(-1,1)
			p2=Rnd(-1,1)
			y2=Rnd(-1,1)
			r2=Rnd(-1,1)
			timer=MilliSecs()
		EndIf
		TurnEntity cube,p,y,r

		l3d = create_3D_line(0,EntityX(cube,1),EntityY(cube,1),EntityZ(cube,1),EntityX(enemy,1),EntityY(enemy,1),EntityZ(enemy,1),255,0,0)
		
		TFormPoint EntityX(enemy,1),EntityY(enemy,1),EntityZ(enemy,1),0,cube
		dx# = TFormedX()
		dy# = TFormedY()
		dz# = TFormedZ()
		rat1# = xsize/Abs(dx)
		rat2# = ysize/Abs(dy)
		rat3# = zsize/Abs(dz)
		rat# = lowest_val(rat1,rat2,rat3)
		nx# = rat * dx
		ny# = rat * dy
		nz# = rat * dz
		TFormPoint nx,ny,nz,cube,0
		PositionEntity hit,TFormedX(),TFormedY(),TFormedZ(),1
		
		PositionEntity spots(free_spot),nx,ny,nz
		free_spot=free_spot+1
		If free_spot>num_spots Then free_spot=1

				

		CameraClsMode cam,1,1
		HideEntity cube
		HideEntity hit
		HideEntity enemy
		ShowEntity l3d
		WireFrame 1
		RenderWorld

		CameraClsMode cam,0,0
		ShowEntity cube
		ShowEntity hit
		ShowEntity enemy
		HideEntity l3d
		WireFrame 0
		RenderWorld
		
		
		frame_time = MilliSecs() - frame_start	
		show_info()

		WaitTimer(fps_timer)
		Flip(1)
		FreeEntity l3d
		If slowmo Then Delay 200
	Wend

	ClearWorld	

	End

Function lowest_val#(v1#,v2#,v3#)

	If v1<v2 And v1<v3
		Return v1
	ElseIf v2<v3
		Return v2
	Else
		Return v3
	EndIf

End Function

;
; Display debug info
;
Function show_info()
	
	If fps_timeout
		frame_count = frame_count + 1

		If MilliSecs() > fps_timeout Then
			fps_timeout = MilliSecs() + 1000 
			fps = frame_count 
			frame_count = 0 
		
			If fps < slowest_fps Or slowest_fps = 0 Then slowest_fps = fps
		EndIf 
		
		If frame_time > slowest_frame Then slowest_frame = frame_time
		
		Color 0,255,0
		Text 10,10," Triangles: " + TrisRendered()
		Color 255,255,0
		Text 10,25," Millisecs: " + frame_time
		Text 10,40,"   Slowest: " + slowest_frame
		Color 0,255,255
		Text 10,55,"       FPS: " + fps
		Text 10,70,"     Worst: " + slowest_fps
		Color 255,255,255
	Else
		; First call initialization.
		fps_timeout = MilliSecs() + 1000 
	EndIf
	
End Function
	
	
Function create_3D_line(mesh,x0#,y0#,z0#,x1#,y1#,z1#,r%=255,g%=255,b%=255) 

	If mesh = 0 
		mesh = CreateMesh() 
		surf = CreateSurface(mesh) 
		EntityFX mesh,1+2+16
	Else 
		last_surf = CountSurfaces(mesh)
		surf = GetSurface(mesh,last_surf)
		If CountVertices(surf) > 30000 Then surf = CreateSurface(mesh)
	End If 

	v0 = AddVertex(surf,x0,y0,z0) 
	v1 = AddVertex(surf,x1,y1,z1)  
	v2 = AddVertex(surf,x0,y0,z0)  
	AddTriangle surf,v0,v1,v2
	
	VertexColor surf,v0,r,g,b
	VertexColor surf,v1,r,g,b
	VertexColor surf,v2,r,g,b

	Return mesh 

End Function


thank you for this demonstration but i can not use tform-points because this is all in 2D. also i need to be able to pick from any point to any direction, and not to the centre. but thanks for your attention and your help attempt :)

No problem. I didn't think it would be much use, but I thought I'd post it just in case.

Sorry, I know nothing about ray tracers. :)


Finding the closest point on a line to another point.
-----------------------------------------------------
(x1, y1, z1), (x2, y2, z2) = ends of line segment
(x3, y3, y3) = point


u =

(x3 - x1)(x2 - x1) + (y3 - y1)(y2 - y1) + (z3 - z1)(z2 - z1)
-----------------------------------------------------------
(x2 - x1)(x2 - x1) + (y2 - y1)(y2 - y1) + (z2 - z1)(z2 - z1)



u = the location between the two end points of the line which is closest to the point.

If u is not between 0 and 1 then the closest point is not between P1 and P2




Find that point on the line closest to the center of the sphere, then calculate the squared distance of it from the center of the sphere. If the squared distance of that is less than the sqaured radius of the sphere, then the line intersects the sphere.


If you need to know the first point where the ray intersects the sphere though, I can't help you there. There's probably some way to calculate that based on the closest point, and the distance of that closest point to the center of the sphere, because the entry and exit points are probably equal distances from that closest point on the line, and those distances are probably related to how deeply the line is embedded in the sphere.

sswift: i already have a sphere intersection routine, but i still need a cube intersection routine. lately i saw a plane routine which i will use soon.

but now i barely understand how sphere picking works!!
thank you for explination :)


edit:

implemented sphere shading. it works simple. the brightness multipler is the vector dot product between the surface normal and a normalised vector that points toward a light!