shadows.mws
Shadows and Projections of 3-D Objects
Sylvain Muise
<smuise@student.math.uwaterloo.ca>
Introduction
This Maple worksheet provides examples and animations of 3-D objects and their shadows. In the first example, a box with two open sides rotates around a light source, and its shadow is projected on the xy plane. In the second example, a rotating ring flies over a hilly landscape, and a travelling sun projects the ring's shadow on the landscape.
Rotating Box
Introduction
This example uses rotation and projection matrices to animate a rotating box and its shadow. The light source is on the z-axis, a little higher than the top of the box. The box rotates around the z-axis, and at the same time rotates in two directions around itself. The shadow of the box is projected onto the xy plane, and the whole is animated for one whole cycle of the box. The animation is presented at the end of the worksheet.
Code
>
restart: with(plots): with(plottools): with(LinearAlgebra):
Warning, the name changecoords has been redefined
Warning, the name arrow has been redefined
>
setoptions3d(scaling=constrained,axes=normal,view=[-23..23,-23..23,0..7]);
The vertices and center of the box are defined:
>
v := [[-2,2,2],[-2,6,2],[2,6,2],[2,2,2],[-2,2,6],[-2,6,6],[2,6,6],[2,2,6]]:
c := [0,4,4]:
The rotation matrices are defined:
The first rotates about the z-axis, the second rotates about the y-axis, and the third rotates about the x-axis.
>
rotMat1 := t -> Matrix([[cos(t),-sin(t),0],[sin(t),cos(t),0],[0,0,1]]):
rotMat2 := t -> Matrix([[cos(t),0,-sin(t)],[0,1,0],[sin(t),0,cos(t)]]):
rotMat3 := t -> Matrix([[1,0,0],[0,cos(t),-sin(t)],[0,sin(t),cos(t)]]):
Let's show how these rotation matrices work. Let's define three vectors and use these matrices on them:
>
setoptions3d(scaling=constrained,view=[-5..5,-5..5,0..6]):
vector1 := Vector([3,2,5]);
vector2 := Vector([4,1,5]);
vector3 := Vector([2,1.5,5.5]);
We can see the Vectors with the plots[arrow] command:
>
arrow1 := plots[arrow](vector1,color=red):
arrow2 := plots[arrow](vector2,color=blue):
arrow3 := plots[arrow](vector3,color=green):
display(arrow1,arrow2,arrow3);
Lets rotate around the z-axis by 3Pi/4 radians by using the first rotation matrix:
>
vector4 := rotMat1(3*Pi/4).vector1:
vector5 := rotMat1(3*Pi/4).vector2:
vector6 := rotMat1(3*Pi/4).vector3:
arrow4 := plots[arrow](vector4,color=orange):
arrow5 := plots[arrow](vector5,color=magenta):
arrow6 := plots[arrow](vector6,color=yellow):
display(arrow1,arrow2,arrow3,arrow4, arrow5, arrow6);
Let's now rotate the vectors around the y and x axes, using rotMat2 and rotMat3. But instead of rotating them around the actual axis as in the previous example, lets make believe the yellow vector is the origin. In order to rotate around vector 6, the yellow one, we have to first translate each of the vectors by -vector6, so that vector 6 is at the origin and the other two are still in the same relative position around vector 6. We then apply whatever rotation matrices we want, then translate them back to their original position by +vector6.
>
vector7 := vector4 - vector6:
vector8 := vector5 - vector6:
vector7 := rotMat2(Pi).vector7: vector7 := rotMat3(Pi).vector7:
vector8 := rotMat2(Pi).vector8: vector8 := rotMat3(Pi).vector8:
vector7 := vector7 + vector6:
vector8 := vector8 + vector6:
arrow7 := plots[arrow](vector7,color=maroon):
arrow8 := plots[arrow](vector8,color=cyan):
display(arrow1,arrow2,arrow3,arrow4,arrow5,arrow6,arrow7,arrow8);
Now the projection part. For simplicity's sake, I have chosen the light source to be on the z-axis. Consider the two following diagrams:
>
line1 := line([0,0,9],[0,27/4,0]): line2 := line([0,0,9],[45/4,27/4,0]): line3 := line([0,0,9],[45/4,0,0]): poly1 := polygon([[0,0,5],[0,3,5],[5,3,5],[5,0,5]],color=red): poly2 := polygon([[0,0,0],[0,27/4,0],[45/4,27/4,0],[45/4,0,0]],color=black):
line4 := line([0,0],[0,9]): line5 := line([0,9],[45/4,0]): line6 := line([0,0],[45/4,0]): line7 := line([0,5],[5,5]): text1 := textplot([-3,5,`d`]): arr1 := arrow([-3,6],[-3,9],.2, .4, .1): arr2 := arrow([-3,4],[-3,0],.2, .4, .1): text2 := textplot([-1.5,7,`d-z`]): text3 := textplot([-1.5,2.5,`z`]): arr3 := arrow([-1.5,7.5],[-1.5,9],.2,.4,.1): arr4 := arrow([-1.5,6.5],[-1.5,5],.2,.4,.1): arr5 := arrow([-1.5,3],[-1.5,5],.2,.4,.1): arr6 := arrow([-1.5,2],[-1.5,0],.2,.4,.1): text4 := textplot([2.5,4,`x`]): text5 := textplot([5.625, -1, `NewX`]): arr7 := arrow([3,4],[5,4],.2,.4,.1): arr8 := arrow([2,4],[0,4],.2,.4,.1): arr9 := arrow([6.25,-1],[11.25,-1],.2,.4,.1): arr10 := arrow([5,-1],[0,-1],.2,.4,.1):
setoptions3d(view=[0..12,0..7,0..10]): setoptions(view=[-3.5..12,-2..10], axes=none):
display(line1,line2,line3,poly1,poly2, orientation=[77,64]);
display(line4,line5,line6,line7,text1, arr1, arr2, text2, text3, arr3, arr4, arr5, arr6, text4, text5, arr7,arr8, arr9, arr10);
Similiar triangles show that
. Therefore
Similarly,
. Where d is the distance from the xy plane to the light source on the z-axis.
So to project, I simply apply these two formulas for the projected x and y, and display them!
>
setoptions3d(scaling=constrained,view=[-5..10,-5..7,0..6]):
vector9 := Vector([vector1[1]/(1-vector1[3]/9), vector1[2]/(1-vector1[3]/9), 0]):
vector10 := Vector([vector2[1]/(1-vector2[3]/9), vector2[2]/(1-vector2[3]/9), 0]):
vector11 := Vector([vector3[1]/(1-vector3[3]/9), vector3[2]/(1-vector3[3]/9), 0]):
arrow9 := plots[arrow](vector9, color=black):
arrow10 := plots[arrow](vector10,color=black):
arrow11 := plots[arrow](vector11,color=black):
display(arrow1, arrow2, arrow3, arrow9, arrow10, arrow11);
Now I apply these procedures to the 8 vertices of the cube previously defined, and animate it!
>
N := 60:
for j from 0 to N do
pVects := map(Vector,v):
pVects := map2(MatrixVectorMultiply,rotMat1(j/N * 2*Pi),pVects):
cVect := Vector(c):
wVect := MatrixVectorMultiply(rotMat1(j/N*2*Pi),cVect):
for i from 1 to 8 do
pVects[i] := rotMat2(j/N*6*Pi) . (pVects[i] - wVect) + wVect:
pVects[i] := rotMat3(j/N*4*Pi) . (pVects[i] - wVect) + wVect:
end do:
p := map(convert,pVects,list):
s := [polygon([p[1], p[2], p[3], p[4]], color=red), polygon([p[1], p[2], p[6], p[5]], color=red), polygon([p[5], p[6], p[7], p[8]], color=red), polygon([p[8], p[7], p[3], p[4]], color=red)]:
sides := display([seq(s[i],i=1..4)]):
cube[j] := [sides]:
for i from 1 to 8 do
shadowPts[i] := [p[i][1]/ (1 - p[i][3]/9),p[i][2]/ (1 - p[i][3]/9),0]:
end do:
shadowS := [polygon([shadowPts[1], shadowPts[2], shadowPts[3], shadowPts[4]],color=black), polygon([shadowPts[1], shadowPts[2], shadowPts[6], shadowPts[5]],color=black), polygon([shadowPts[5], shadowPts[6], shadowPts[7], shadowPts[8]],color=black), polygon([shadowPts[8], shadowPts[7], shadowPts[3], shadowPts[4]],color=black)]:
shadowSides := display([seq(shadowS[i],i=1..4)]):
shadow[j] := [shadowSides]:
end do:
>
ani1 := display(seq(cube[i],i=0..N),insequence=true):
>
ani2 := display(seq(shadow[i],i=0..N),insequence=true):
>
setoptions3d(scaling=constrained,axes=normal,view=[-20..20,-20..20,0..9]);
Animation
>
display([ani1,ani2]);
>
Rotating Ring
Introduction
In this example, a floating and rotating ring travels over a hilly landscape, while the sun traverses the sky in its usual journey. The resulting shadow from the sun and ring is shown on the landscape.
Code
>
restart: with(plots): with(plottools): with(LinearAlgebra):
Warning, the name changecoords has been redefined
Warning, the name arrow has been redefined
>
sun := T -> [10*cos(T),3.5,10*sin(T)]:
>
rotMat := u -> Matrix([[cos(u),0,-sin(u)],[0,1,0],[sin(u),0,cos(u)]]):
>
w := u -> Vector([5*u + 5,-5*u,2]):
>
ringRot := (t,u) -> Vector([cos(t)+5*u+5, sin(t)-5*u,2]):
>
ring := (t,u) -> convert(rotMat(u) . (ringRot(t,u/(-2*Pi)) - w(u/(-2*Pi))) + w(u/(-2*Pi)),list):
>
f := (x,y) -> sin(x)*sin(y):
>
plot1 := plot3d(f(x,y),x=-3..8,y=-5..9, style=patchnogrid):
>
ray := [x,y,z] - sun(T) = expand(s*( sun(T) - ring(t,u))):
>
ray := [solve(lhs(ray)[1] = rhs(ray)[1],x), solve(lhs(ray)[2] = rhs(ray)[2],y), solve(lhs(ray)[3] = rhs(ray)[3],z)]:
>
N := 30:
>
ringPlot := j -> spacecurve(ring(t,(j/N) * (2*Pi)),t=0..2*Pi,color=red, thickness=3):
>
ringMove := display(seq(ringPlot(j),j=0..N),insequence=true):
>
sunPlot := j -> sphere(sun(Pi/6 + ((j/N)*2*Pi/3)), 2, style=patchnogrid, color=yellow):
>
sunMove := display(seq(sunPlot(j),j=0..N),insequence=true):
>
rayLine := subs(T = Pi/6,ray):
rayLine := subs(u = 0, rayLine):
n := 15:
for i from 1 to n do
rays := subs(t = i * 2 * Pi/n, rayLine):
lines[i] := line(sun(Pi/6), subs(s=fsolve(f(rays[1],rays[2]) = rays[3], s),rays)):
end do:
>
for j from 0 to N do
temp2 := subs(T = Pi/6 + ((j/N) * 2*Pi/3), ray):
temp3 := subs(u = (j/N) * (2*Pi), temp2):
n := 30:
for i from 1 to n do
temp := subs(t = i * 2 * Pi/n, temp3):
points[i] := subs(s=fsolve(f(temp[1],temp[2]) = temp[3], s),temp):
end do:
shadow[j] := curve([seq(points[i],i=1..n)],color=black,thickness=3):
end do:
>
ani := display(seq(shadow[j],j=0..N), insequence=true):
Graphic and Animation
To draw the shadow of the ring, lines are computed from the sun to a number of points on the ring, and then they are extended until they meet the surface, as shown in the following image. (In the image, only 15 lines are drawn for clarity, but in reality 30 are computed.) The shadow is then drawn by joining all these points, and the animation is created by doing this for every frame of time, in which both the sun and ring move.
>
display([sunPlot(0),ringPlot(0),plot1, seq(lines[i],i=1..n)], view=[0..11, -5..6, -2..7], scaling=constrained);
Here is the complete animation:
>
display([plot1,ani,sunMove,ringMove],scaling=constrained,view=[-8..10,-5..9,-2..10]);
>