[bb] ray->plane ray->triangle intersection by elias_t [ 1+ years ago ]

[BlitzBot]

Title : ray->plane ray->triangle intersection
Author : elias_t
Posted : 1+ years ago

Description : ;2 functions which find out if a line segment and a plane or triangle intersect
;and return the intersection point.
The ray->tri function is a mod of sswift's posted code.

Code :

Code: BlitzBasic

1. ;FIRST FILE.
2.  
3. ;-------------------------------------------------------------------------
4. ;A function that finds out if a line segment and a plane intersect
5. ;and calculates the intersection point in the array picked(2)
6. ;where picked(0)=pickedx,..1=y,..2=z
7. ;
8. ;After a tutorial by Paul Bourke.
9. ;
10. ;by elias_t
11. ;
12. ;-------------------------------------------------------------------------
13.  
14. ;stores the picked location
15. Dim picked#(2)
16.  
17.  
18.  
19. Function ray_plane(p1x#,p1y#,p1z#, p2x#,p2y#,p2z#, pax#,pay#,paz#, pbx#,pby#,pbz#, pcx#,pcy#,pcz#)
20.  
21. Local d#
22. Local total#,denom#,mu#
23. Local nx#,ny#,nz#
24.  
25.                                
26. Local dx1#=(pcx - pax)
27. Local dy1#=(pcy - pay)
28. Local dz1#=(pcz - paz)
29.  
30. Local dx2#=(pbx - pax)
31. Local dy2#=(pby - pay)
32. Local dz2#=(pbz - paz)
33.  
34. nx#=(dy1*dz2)-(dz1*dy2)
35. ny#=(dz1*dx2)-(dx1*dz2)
36. nz#=(dx1*dy2)-(dy1*dx2)
37.  
38. ;------------------------------------
39.  
40.  
41.  
42.    d = - nx * pax - ny * pay - nz * paz
43.  
44.    ;Calculate the position on the Line that intersects the plane
45.    denom = nx * (p2x - p1x) + ny * (p2y - p1y) + nz * (p2z - p1z);
46.  
47.    If (Abs(denom) < 0.0001) Return 0;Line And plane don't intersect
48.      
49.    mu = - (d + nx * p1x + ny * p1y + nz * p1z) / denom
50.    picked(0) = p1x + mu * (p2x - p1x)
51.    picked(1) = p1y + mu * (p2y - p1y)
52.    picked(2) = p1z + mu * (p2z - p1z)
53.        
54.         ;comment this out if you want an infinite ray
55.    If (mu < 0 Or mu > 1) Return 0;Intersection Not along Line segment
56.      
57.    Return 1
58.  
59. End Function
60.  
61.  
62.  
63.  
64.  
65. ;-------------------------------------------------------------------------
66. ;EXAMPLE
67. ;quick and messy but works ok.
68.  
69. Graphics3D 640,480,32,2
70.  
71. cam=CreateCamera()
72. CameraRange cam,1,20000
73. PositionEntity cam,0,5,-10
74. li=CreateLight()
75.  
76. ;start and end point of line segment
77. p1x#=0:p1y#=10:p1z#=-10
78. p2x#=-120:p2y#=-500:p2z#=1000
79.  
80.  
81.  
82. ;visualize the triangle that defines a plane
83. m=CreateMesh()
84. s=CreateSurface(m)
85. v0=AddVertex(s,1,0,0)
86. v1=AddVertex(s,0,0,1)
87. v2=AddVertex(s,-1,0,0)
88. v3=AddVertex(s,0,0,-1);extra vert to make a plane
89. AddTriangle (s,v0,v2,v1)
90.  
91. EntityColor m,0,255,0
92. UpdateNormals m
93.  
94.  
95.  
96.  
97. RotateMesh m,20,0,30
98. ScaleMesh m,400,400,400
99.  
100. ;3 vertexes of the triangle that define a plane
101. v0x#=VertexX(s,v0):v0y#=VertexY(s,v0):v0z#=VertexZ(s,v0)
102. v1x#=VertexX(s,v1):v1y#=VertexY(s,v1):v1z#=VertexZ(s,v1)
103. v2x#=VertexX(s,v2):v2y#=VertexY(s,v2):v2z#=VertexZ(s,v2)
104.  
105.  
106.  
107.  
108. ;create the visual ray
109. test=CopyMesh(m)
110. EntityColor test,255,255,0
111. EntityFX test,16
112.  
113. s1=GetSurface(test,1)
114.  
115. ;the +-.2 is for visualisation purposes
116. VertexCoords(s1,0,p1x+.2,p1y,p1z)
117. VertexCoords(s1,1,p2x,p2y,p2z)
118. VertexCoords(s1,2,p1x-.2,p1y,p1z)
119.  
120.  
121. PositionEntity cam,0,15,-10
122. PointEntity cam ,m
123.  
124.  
125. ;intersction pointer
126. cc=CreateCube()
127. EntityColor cc,255,0,0
128. ScaleEntity cc,.2,.2,.2
129.  
130. ;add the extra tri to make a plane
131. AddTriangle (s,v0,v3,v2)
132. UpdateNormals m
133.  
134. While Not KeyHit(1)
135.  
136. RenderWorld()
137.  
138. If KeyDown(200) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1)+1,VertexZ(s1,1))
139. If KeyDown(208) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1)-1,VertexZ(s1,1))
140.  
141. If KeyDown(203) Then VertexCoords(s1,1,VertexX(s1,1)-1,VertexY(s1,1),VertexZ(s1,1))
142. If KeyDown(205) Then VertexCoords(s1,1,VertexX(s1,1)+1,VertexY(s1,1),VertexZ(s1,1))
143.  
144. If KeyDown(201) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1),VertexZ(s1,1)+1)
145. If KeyDown(209) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1),VertexZ(s1,1)-1)
146.  
147.  
148. dx#=VertexX(s1,1)
149. dy#=VertexY(s1,1)
150. dz#=VertexZ(s1,1)
151.  
152. res=ray_plane(p1x,p1y,p1z, dx,dy,dz  ,v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z)
153. Text 10,10,res
154.  
155.  
156. If res
157. PositionEntity cc,picked(0),picked(1),picked(2),1
158. Text 10,100,picked(0)
159. Text 10,120,picked(1)
160. Text 10,140,picked(2)
161. EndIf
162.  
163.  
164.  
165. Flip
166.  
167. Wend
168.  
169.  
170. End
171.  
172.  
173.  
174.  
175.  
176. ;SECOND FILE--------------------------------------------------------
177.  
178.  
179. ;-------------------------------------------------------------------------
180. ;A function that finds out if a line segment and a triangle intersect
181. ;and calculates the intersection point in the array picked(2)
182. ;where picked(0)=pickedx,..1=y,..2=z
183. ;
184. ;This actually the function posted by sswift and originally made by Tomas Moller.
185. ;Only here we get also the intersection point !
186. ;
187. ;mod by elias_t
188. ;
189. ;-------------------------------------------------------------------------
190.  
191. ;stores the picked location
192. Dim picked#(2)
193.  
194.  
195. Function rit(Px#,Py#,Pz#, Dx#,Dy#,Dz#, V0x#,V0y#,V0z#, V1x#,V1y#,V1z#, V2x#,V2y#,V2z#, Extend_To_Infinity=1, Cull_Backfaces=0)
196.        
197.         ;-
198.         E1x# = V2x# - V0x#
199.         E1y# = V2y# - V0y#
200.         E1z# = V2z# - V0z#
201.  
202.         ;-
203.         E2x# = V1x# - V0x#
204.         E2y# = V1y# - V0y#
205.         E2z# = V1z# - V0z#
206.  
207.         ; Hxyz = Crossproduct(Dxyz, E2xyz)
208.         Hx# = (Dy# * E2z#) - (E2y# * Dz#)
209.         Hy# = (Dz# * E2x#) - (E2z# * Dx#)
210.         Hz# = (Dx# * E2y#) - (E2x# * Dy#)
211.  
212.         ; Calculate the dot product of the above vector and the vector between point 0 and point 2.
213.         A# = (E1x# * Hx#) + (E1y# * Hy#) + (E1z# * Hz#)
214.  
215.         ;cull
216.         If (Cull_Backfaces = 1) And (A# >= 0) Then Return 0
217.                
218.         ;parralel
219.         If (A# > -0.00001) And (A# < 0.00001) Then Return 0
220.        
221.         ;Inverse Determinant
222.         F# = 1.0 / A#
223.  
224.         ;-
225.         Sx# = Px# - V0x#
226.         Sy# = Py# - V0y#
227.         Sz# = Pz# - V0z#
228.        
229.         ; U# = F# * (DotProduct(Sxyz, Hxyz))
230.         U# = F# * ((Sx# * Hx#) + (Sy# * Hy#) + (Sz# * Hz#))
231.        
232.         ;check u
233.         If (U# < 0.0) Or (U# > 1.0) Return 0
234.  
235.         ; Qxyz = CrossProduct(Sxyz, E1xyz)
236.         Qx# = (Sy# * E1z#) - (E1y# * Sz#)
237.         Qy# = (Sz# * E1x#) - (E1z# * Sx#)
238.         Qz# = (Sx# * E1y#) - (E1x# * Sy#)
239.        
240.         ; V# = F# * DotProduct(Dxyz, Qxyz)
241.         V# = F# * ((Dx# * Qx#) + (Dy# * Qy#) + (Dz# * Qz#))
242.  
243.         ;check v
244.         If (V# < 0.0) Or ((U# + V#) > 1.0) Return 0
245.  
246.  
247.         T# = F#*((E2x#*Qx#)+(E2y#*Qy#)+(E2z#*Qz#))
248.        
249.         If T#<0 Return 0
250.        
251.  
252.         If Extend_To_Infinity=0 And T#>1 Return 0
253.  
254.  
255.  
256.  
257. ;-------------------------------------------------
258. ;Calculate intersection point
259. nx#=(E1y*E2z)-(E1z*E2y)
260. ny#=(E1z*E2x)-(E1x*E2z)
261. nz#=(E1x*E2y)-(E1y*E2x)
262. d# = -  nx*V0x - ny*V0y - nz*V0z
263.  
264. denom# = nx*Dx + ny*Dy + nz*Dz
265. mu# = - (d + nx*Px + ny*Py + nz*Pz) / denom
266.  
267. picked(0) = Px + mu * DX
268. picked(1) = Py + mu * Dy
269. picked(2) = Pz + mu * Dz
270. ;-------------------------------------------------
271.  
272.  
273.  
274.         ;intersects            
275.         Return 1
276.  
277. End Function
278.  
279.  
280.  
281.  
282.  
283.  
284.  
285.  
286. ;Example
287. Graphics3D 640,480,32,2
288.  
289. cam=CreateCamera()
290. PositionEntity cam,0,5,-10
291. li=CreateLight()
292.  
293. ;start and end point of line segment
294. p1x#=0:p1y#=0:p1z#=-10
295. p2x#=0:p2y#=5:p2z#=10
296.  
297. ;3 vertexes of the triangle
298. v0x#=8:v0y#=0:v0z#=-1
299. v1x#=0:v1y#=5:v1z#=0
300. v2x#=-5:v2y#=0:v2z#=1
301.  
302.  
303. dx#=p2x-p1x
304. dy#=p2y-p1y
305. dz#=p2z-p1z
306.  
307. st=MilliSecs()
308. For i=1 To 1000000
309. res=rit(p1x,p1y,p1z,dx,dy,dz,v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z)
310. Next
311. et=MilliSecs()
312. dt=(et-st)
313.  
314.  
315. m=CreateMesh()
316. s=CreateSurface(m)
317. v0=AddVertex(s,v0x,v0y,v0z)
318. v1=AddVertex(s,v1x,v1y,v1z)
319. v2=AddVertex(s,v2x,v2y,v2z)
320. AddTriangle (s,v0,v2,v1)
321. EntityColor m,0,255,0
322. UpdateNormals m
323.  
324. test=CopyMesh(m)
325. EntityColor test,255,255,0
326. EntityFX test,16
327.  
328. s1=GetSurface(test,1)
329.  
330. ;the +-.2 is for visualisation purposes
331. VertexCoords(s1,0,p1x+.2,p1y,p1z)
332. VertexCoords(s1,1,p2x,p2y,p2z)
333. VertexCoords(s1,2,p1x-.2,p1y,p1z)
334.  
335.  
336. PositionEntity cam,0,15,-10
337. PointEntity cam ,m
338.  
339.  
340. ;intersction pointer
341. cc=CreateCube()
342. EntityColor cc,255,0,0
343. ScaleEntity cc,.2,.2,.2
344.  
345.  
346.  
347.  
348.  
349. While Not KeyHit(1)
350.  
351. RenderWorld()
352.  
353. If KeyDown(200) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1)+.1,VertexZ(s1,1))
354. If KeyDown(208) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1)-.1,VertexZ(s1,1))
355.  
356. If KeyDown(203) Then VertexCoords(s1,1,VertexX(s1,1)-.1,VertexY(s1,1),VertexZ(s1,1))
357. If KeyDown(205) Then VertexCoords(s1,1,VertexX(s1,1)+.1,VertexY(s1,1),VertexZ(s1,1))
358.  
359. If KeyDown(201) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1),VertexZ(s1,1)+.1)
360. If KeyDown(209) Then VertexCoords(s1,1,VertexX(s1,1),VertexY(s1,1),VertexZ(s1,1)-.1)
361.  
362.  
363. dx#=VertexX(s1,1)-p1x
364. dy#=VertexY(s1,1)-p1y
365. dz#=VertexZ(s1,1)-p1z
366.  
367. res=rit(p1x,p1y,p1z,dx,dy,dz,v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,0)
368.  
369. Text 10,10,res
370.  
371. If res=1
372. PositionEntity cc,picked(0),picked(1),picked(2),1
373. Text 10,100,picked(0)
374. Text 10,120,picked(1)
375. Text 10,140,picked(2)
376. EndIf
377.  
378. Text 10,40,dt
379.  
380. Flip
381.  
382. Wend
383.  
384.  
385. End

Comments : none...