This code calculates and displays the intersection lines of three self-intersecting tori:
sync on
type vec3
x as float
y as float
z as float
endtype
global vec3_result as vec3
n as vec3
p as vec3
n.x = 1
n.y = 2
n.z = 3
vec3_setlen(n, 1.0)
n = vec3_result
autocam off
zero as vec3
p1 as vec3
p2 as vec3
p3 as vec3
n12 as vec3
n23 as vec3
n31 as vec3
c12 as vec3
c23 as vec3
c31 as vec3
p2.x = 10.0
p3.z = 10.0
a12# = 20.0
a23# = 25.0
a31# = 30.0
vec3_add(p1, p2)
vec3_mul(vec3_result, 0.5)
c12 = vec3_result
vec3_sub(p2, p1)
l# = vec3_setlen(vec3_result, 1.0)
n12 = vec3_result
r1_12# = inscribed_r1(l#, a12#)
r2_12# = inscribed_r2(l#, a12#)
vec3_add(p2, p3)
vec3_mul(vec3_result, 0.5)
c23 = vec3_result
vec3_sub(p3, p2)
l# = vec3_setlen(vec3_result, 1.0)
n23 = vec3_result
r1_23# = inscribed_r1(l#, a23#)
r2_23# = inscribed_r2(l#, a23#)
vec3_add(p3, p1)
vec3_mul(vec3_result, 0.5)
c31 = vec3_result
vec3_sub(p1, p3)
l# = vec3_setlen(vec3_result, 1.0)
n31 = vec3_result
r1_31# = inscribed_r1(l#, a31#)
r2_31# = inscribed_r2(l#, a31#)
make object box 1, 0.5, 0.5, 0.5
make object box 2, 0.5, 0.5, 0.5
make object box 3, 0.5, 0.5, 0.5
hide object 1
hide object 2
hide object 3
set object emissive 1, 0xFFFF0000
set object emissive 2, 0xFF00FF00
set object emissive 3, 0xFFFFFF00
set object diffuse 1, 0
set object diffuse 2, 0
set object diffuse 3, 0
do
position camera 0, 0, 0
rotate camera xang#, yang#, 0.0
move camera -100.0
inc xang#, mousemovey()*0.5
inc yang#, mousemovex()*0.5
for i = 1 to 5
p.x = random(100.0)-50.0
p.y = random(100.0)-50.0
p.z = random(100.0)-50.0
vec3_result = p
for j = 1 to 100
vec3_project_onto_torus(vec3_result, c12, n12, r1_12#, r2_12#)
vec3_project_onto_torus(vec3_result, c23, n23, r1_23#, r2_23#)
next j
p = vec3_result
id = find free object()
instance object id, 1
position object id, p.x, p.y, p.z
next i
for i = 1 to 5
p.x = random(100.0)-50.0
p.y = random(100.0)-50.0
p.z = random(100.0)-50.0
vec3_result = p
for j = 1 to 100
vec3_project_onto_torus(vec3_result, c23, n23, r1_23#, r2_23#)
vec3_project_onto_torus(vec3_result, c31, n31, r1_31#, r2_31#)
next j
p = vec3_result
id = find free object()
instance object id, 2
position object id, p.x, p.y, p.z
next i
for i = 1 to 5
p.x = random(100.0)-50.0
p.y = random(100.0)-50.0
p.z = random(100.0)-50.0
vec3_result = p
for j = 1 to 100
vec3_project_onto_torus(vec3_result, c31, n31, r1_31#, r2_31#)
vec3_project_onto_torus(vec3_result, c12, n12, r1_12#, r2_12#)
next j
p = vec3_result
id = find free object()
instance object id, 3
position object id, p.x, p.y, p.z
next i
sync
loop
function inscribed_r1(l as float, a as float)
v# = l*tan(90.0-a)
endfunction v#
function inscribed_r2(l as float, a as float)
v# = l/cos(90.0-a)
endfunction v#
function vec3_project_onto_torus(a as vec3, c as vec3, n as vec3, r1 as float, r2 as float)
vec3_project_onto_circle(a, c, n, r1)
vec3_project_onto_sphere(a, vec3_result, r2)
endfunction
function vec3_project_onto_circle(a as vec3, c as vec3, n as vec3, r as float)
vec3_project_onto_plane(a, c, n)
vec3_project_onto_sphere(vec3_result, c, r)
endfunction
function vec3_project_onto_plane(a as vec3, p as vec3, n as vec3)
vec3_sub(p, a)
vec3_mul(n, vec3_dot(vec3_result, n))
vec3_add(vec3_result, a)
endfunction
function vec3_project_onto_sphere(a as vec3, c as vec3, r as float)
vec3_sub(a, c)
vec3_setlen(vec3_result, r)
vec3_add(vec3_result, c)
endfunction
function vec3_setlen(a as vec3, l as float)
ol# = vec3_length(a)
vec3_mul(a, l/ol#)
endfunction ol#
function vec3_length(a as vec3)
v# = sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
endfunction v#
function vec3_dot(a as vec3, b as vec3)
v# = a.x*b.x + a.y*b.y + a.z*b.z
endfunction v#
function vec3_mul(a as vec3, b as float)
vec3_result.x = a.x*b
vec3_result.y = a.y*b
vec3_result.z = a.z*b
endfunction
function vec3_add(a as vec3, b as vec3)
vec3_result.x = a.x+b.x
vec3_result.y = a.y+b.y
vec3_result.z = a.z+b.z
endfunction
function vec3_sub(a as vec3, b as vec3)
vec3_result.x = a.x-b.x
vec3_result.y = a.y-b.y
vec3_result.z = a.z-b.z
endfunction
As well as producing some nice colourful curves, the results are actually very useful. I originally came to write this because I was trying to calculate the position in space of a video camera based on three known points in its field of view.
This problem turns out to be very complicated, but after much doing of maths I managed to work out that it was equivalent to the points of intersection of three tori.
Unfortunately calculating the points where three tori intersect is also very difficult (requires solving three simultaneous non-linear equations) but I eventually found this iterative numerical solution which gives very good results (as you can see).
Also, the vector functions there are very useful on their own.
Maths stuff:
The original problem was to use three points whose 2d positions are known on the image produced by a camera to work out the location of that camera.
The first problem is that the locations of the points are not only affected by the camera's position, but also by its rotating. The rotation and position each have 3 components giving a total 6 unknowns to be solved simultaneously based on 6 input values (3 X,Y pairs), which is not very practical.
To get around this I found a way to separate the position and rotation information, simplifying the problem to 3 unknowns: work out the angle between each pair of points (this is approximately proportional to the distance between each pair, but can be worked out exactly without much difficulty).
However the camera is rotated, the angle between each pair of points will not change. Now the problem is how to convert the 3 angles into a position.
In 2d, the set of points which make a fixed angle with two other points always make a circle:
inscribed angle theorem. The special case of this which is most well known is that if you pick any point on a semicircle and draw lines to the end points of that semicircle, the angle between them is always 90 degress.
It turns out to be very simple to extend this to 3d: instead of a circle being produced, a circle rotated through another circle is produced (a torus). Each pair of points and the angle between them produces a torus, so now the problem is to find a point that exists on all 3 tori (ie. their intersection).
There's no simple way to do this, but it turns out to be fairly easy to project a random point in space onto the nearest point on a torus, (see the function "vec3_project_onto_torus"), so if points are picked at random and then repeatedly projected onto each torus, they will very quickly end up at the points where they intersect.
Displaying a number of discrete points is not particularly interesting though, so the program instead draws the intersection between each pair of tori in different colours, and the points where those three intersect are the solutions.
[b]
