Application Center - Maplesoft

# Shadows and projections of 3-D objects

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

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:

end do:

> ani1 := display(seq(cube[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:
end do: