Application Center - Maplesoft

App Preview:

Shadows and projections of 3-D objects

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

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]);

vector1 := _rtable[11830676]

vector2 := _rtable[11824920]

vector3 := _rtable[11816044]

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);

[Maple Plot]

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);

[Maple Plot]

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);

[Maple Plot]

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);

[Maple Plot]

[Maple Plot]

Similiar triangles show that NewX/d = x/(d-z) . Therefore NewX = x/(1-z/d)

Similarly, NewY = y/(1-z/d) . 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);

[Maple Plot]

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]);

[Maple Plot]

>

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);

[Maple Plot]

Here is the complete animation:

> display([plot1,ani,sunMove,ringMove],scaling=constrained,view=[-8..10,-5..9,-2..10]);

[Maple Plot]

>