Putting a Scene in Three Dimensions onto a Two-Dimensional Canvas

Introduction

OBJECTIVE: Learn to use linear transformations to map points from three dimensions onto a two-dimensional space.

Since time immemorial, artists have faced a problem that computer graphics artists are facing today, and that is how to picture three-dimensional objects on a two-dimensional plane. The representation of the objects may vary, depending on the position of the "eye of the beholder." In the project that follows, we analyze two aspects of this problem. Both involve mapping points in three-dimensional space to a plane. Part I addresses a problem in which the viewpoint belongs to a single beholder. The remainder of the lab focuses on parallel projections that resemble a situation in which the results of X-rays demonstrate the need for a CAT-Scan. The linear algebra behind such mappings is introduced in Part V.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: Mapping to the Y-Z Plane from a Single Point

First, consider the problem from the text, where all points are mapped to the y-z plane. The viewer is at the point ((Subscript[x, 0]), 0, 0), and any point (x, y, z) between the viewer and the y-z plane appears at the point (0, t y, t z) on the y-z plane, where t will be a function of the viewer's position and the x-coordinate of the point in 3-space. Verify that t = ((Subscript[x, 0])/((Subscript[x, 0]) - x)).
The following commands demonstrate this mapping by computing image points as specified after randomly generating points in the domain that lie between the observer and the y-z plane.

Clear[x,y,z,x0,t,newx,newy,newz,domain,image]
Print["The number of domain and image points is ",numberofpoints=1000]
x0=500;
t[x_,x0_]:=(x0/(x0-x))
newx[x_,y_,z_]:=0
newy[x_,y_,z_]:=t[x,x0] y
newz[x_,y_,z_]:=t[x,x0] z
domain=Table[{RandomReal[{0,100}],RandomReal[{-100,100}],RandomReal[{-100,100}]},{i,1,numberofpoints}];
image=Table[{newx[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newy[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newz[domain[[i,1]],domain[[i,2]],domain[[i,3]]]},{i,1,Length[domain]}];
partimage=Table[image[[i]],{i,1,10}];
PrependTo[partimage,{"x coordinates","y coordinates","z coordinates"}];
Print["The first ten image points are: ",TableForm[partimage]]

To view the points in the image, we  use the command ListPlot3D.

sp=ListPointPlot3D[image,PlotStyle->PointSize[0.01`]]

Let's suppose that there is an object in the region bounded approximately by z=100-0.5(Sqrt[(x^2)+(y^2)]) and z=0.005((x^2)+(y^2)).  We will first plot the object in three-space and then separate out the image points that have been mapped from that region and plot them in pink.

Clear[r]
p1=RevolutionPlot3D[100-0.5` r,{r,0,100},{θ,0,2 π},DisplayFunction->Identity];
p2=RevolutionPlot3D[0.005` (r^2),{r,0,100},{θ,0,2 π},DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction,ViewPoint->{-0.637`,-0.917`,-0.018`}]
listout={};
listin={};
Do[If[domain[[i,3]]<100-0.5` √(((domain[[i,1]])^2)+((domain[[i,2]])^2))&&domain[[i,3]]>0.005` (((domain[[i,1]])^2)+((domain[[i,2]])^2)),AppendTo[listin,{image[[i,1]],image[[i,2]],image[[i,3]]}],AppendTo[listout,{image[[i,1]],image[[i,2]],image[[i,3]]}]],{i,1,Length[domain]}]
pout=ListPointPlot3D[listout,DisplayFunction->Identity];
pin=ListPointPlot3D[listin,PlotStyle->{RGBColor[1,0,1],PointSize[0.02`]},DisplayFunction->Identity];
Show[pin,pout,DisplayFunction->$DisplayFunction]

You Try It: Part I

If you had a parallel projection onto the y-z plane, any (x, y, z) would map to (0, y, z). In the above problem, let (Subscript[x, 0]) get larger and larger, and compare your image points to your domain points to see if you are approaching this state. Do this by making the number in red bigger and bigger.

Clear[x,y,z,x0,t,newx,newy,newz,domain,image]
Print["The number of domain and image points is ",numberofpoints=100]
x0=10000;t[x_,x0_]:=(x0/(x0-x))
newx[x_,y_,z_]:=0
newy[x_,y_,z_]:=t[x,x0] y
newz[x_,y_,z_]:=t[x,x0] z
domain=Table[{RandomReal[{0,100}],RandomReal[{-100,100}],RandomReal[{-100,100}]},{i,1,numberofpoints}];
image=Table[{newx[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newy[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newz[domain[[i,1]],domain[[i,2]],domain[[i,3]]]},{i,1,Length[domain]}];
Clear[list]
list=Table[{domain[[i,j]],image[[i,j]]},{i,1,10},{j,1,3}];
PrependTo[list,{"x coordinates","y coordinates","z coordinates"}];
Print["The first ten domain and image points are:",TableForm[list]]

Are your image values for y and z close to your domain values for y and z respectively?

Part  II: Parallel Mapping to Any Plane Using a Linear Transformation - Identifying the Image

We now consider another type of linear mapping onto a plane. Instead of projections from a single point, we consider parallel projections onto a plane. You could think of this as a generalization of the above case, where the beholder moves farther and farther away ((Subscript[x, 0])->∞). Your domain is the set of all points in three-dimensional space. Begin by finding the image for several points in your domain, and use any three of these distinct noncollinear points in your image to write the equation of the plane in which they lie. Verify that other image points you found satisfy the equation of this plane, and then plot the image to verify visually that the image is a plane.

The transformation described below is defined as follows: 
newx   =x +   y + 2z
newy  =2x +   y  -   z 	
newz  =  3x + 2y  +  z
Other transformations will work, provided that exactly one of the right-hand functions can be written as a linear combination of the other two. In this case, the first right-hand side function equals the third one minus the second one. This restriction guarantees that all your image points will lie in a plane.
Without loss of generality, we focus on domain points with integer coordinates between -100 and 100. We use integers to avoid making roundoff errors when we check to see that image points lie on the plane. One-thousand image points are computed, and all are plotted, but only the first 10 are listed.

Clear[x,y,z,newx,newy,newz,domain,image]
newx[x_,y_,z_]:=x+y+2 z
newy[x_,y_,z_]:=2 x+y-z
newz[x_,y_,z_]:=3 x+2 y+z
Print["The number of domain and image points is ",numberofpoints=1000]
domain=Table[RandomInteger[{-100,100},3],{i,1,numberofpoints}];
image=Table[{newx[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newy[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newz[domain[[i,1]],domain[[i,2]],domain[[i,3]]]},{i,1,Length[domain]}];
Clear[partimage]
partimage=Table[image[[i]],{i,1,10}];
PrependTo[partimage,{"x coordinates","y coordinates","z coordinates"}];
Print["The first ten image points are: "]
Print[TableForm[partimage]]

To determine the equation of the plane in which the first three image points lie, we specify two vectors in the plane and compute their cross product. Knowing that one point in the image is {0, 0, 0} guarantees that the plane passes through the origin.

vector1=image[[1]]-image[[2]];
vector2=image[[1]]-image[[3]];
normal=Cross[vector1,vector2];
Print["The normal to the plane is ",normaltoplane=normal/normal[[1]]]
Clear[x,y,z]
Print["The equation of the plane onto which the points are projected is ",equationofplane[x_,y_,z_]=normaltoplane[[1]]x+normaltoplane[[2]]y+normaltoplane[[3]]z==0]

Check to see if all points you found in the domain satisfy the equation of the plane. We start our iterator at 0 and then increase it by one each time an image point satisfies the equation of the plane.

works=0;
Do[If[equationofplane[image[[i,1]],image[[i,2]],image[[i,3]]],works=works+1],{i,1,Length[image]}] 
works

It looks as though, as predicted, every point in the image lies on the specified plane.

View the points in the image.

sp=ListPointPlot3D[image,PlotStyle->PointSize[0.01`]]

You Try It: Part II

Try the transformation described below: 
newx   =  2x + 2y +   z 	
newy  =  4x + 4y  - 7z 	
newz  =  3x + 3y  + 9z
Change the expressions in red.

Clear[x,y,z,newx,newy,newz,domain,image]
newx[x_,y_,z_]:=1x+1y+2 z
newy[x_,y_,z_]:=2 x+1y-1z
newz[x_,y_,z_]:=3 x+2 y+1z
Print["The number of domain and image points is ",numberofpoints=1000]
domain=Table[RandomInteger[{-100,100},3],{i,1,numberofpoints}];
image=Table[{newx[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newy[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newz[domain[[i,1]],domain[[i,2]],domain[[i,3]]]},{i,1,Length[domain]}];
Clear[partimage]
partimage=Table[image[[i]],{i,1,10}];
PrependTo[partimage,{"x coordinates","y coordinates","z coordinates"}];
Print["The first ten image points are: "]
Print[TableForm[partimage]]

Find the plane in which the image points lie.

vector1=image[[1]]-image[[2]];
vector2=image[[1]]-image[[3]];
normal=Cross[vector1,vector2];
Print["The normal to the plane is ",normaltoplane=normal/normal[[1]]]
Clear[x,y,z]
Print["The equation of the plane onto which the points are projected is ",equationofplane[x_,y_,z_]=normaltoplane[[1]]x+normaltoplane[[2]]y+normaltoplane[[3]]z==0]

Check to see if all points you found in the domain satisfy the equation of the plane. Start your iterator at 0 and then increase it by one each time an image point satisfies the equation of the plane.

works=0;
Do[If[equationofplane[image[[i,1]],image[[i,2]],image[[i,3]]],works=works+1],{i,1,Length[image]}] 
works

Does every point in the image lie on the plane specified?

View the points in the image.

sp=ListPointPlot3D[image,PlotStyle->PointSize[0.01`]]

Do your points look as though they lie on a plane? If not, can you see where there is a problem?

Part III: Analyzing the Mapping

We return to the first mapping defined in Part II. Identify the points in the domain that map to the origin O{0, 0, 0}, and verify that they all lie along the same line. Write the equations for this line.

Off[Solve::svars]
Clear[x,y,z,newx,newy,newz,domain,image]
newx[x_,y_,z_]:=x+y+2z
newy[x_,y_,z_]:=2x+y-z
newz[x_,y_,z_]:=3x+2y+z
maptoO=Solve[{newx[x,y,z]==0,newy[x,y,z]==0,newz[x,y,z]==0},{x,y,z}]

Now select another point in the image, say P{10,-20,-10}, and find the domain points that map to P. Note that they all fall along a line, and observe how this line is related to the line of points that map to the origin.

maptoP=Solve[{newx[x,y,z]==10,newy[x,y,z]==-20,newz[x,y,z]==-10},{x,y,z}]

Check this out for another point in the image Q(-40,70,30) and see if you can make a generalization about the mapping in terms of sets of points in the domain that map to a given point in the image.

maptoQ=Solve[{newx[x,y,z]==-40,newy[x,y,z]==70,newz[x,y,z]==30},{x,y,z}]

Part IV: Application to Medical X-Ray and CAT-Scan Technologies

We will now demonstrate how points within a tumor would show on an X-ray slide.

Clear[x,y,z,newx,newy,newz,domain,image,listin,listout,pin,pout]
newx[x_,y_,z_]:=x+y+2 z
newy[x_,y_,z_]:=2 x+y-z
newz[x_,y_,z_]:=3 x+2 y+z
Print["The number of domain and image points is ",numberofpoints=1000]
domain=Table[RandomInteger[{-100,100},3],{i,1,numberofpoints}];
image=Table[{newx[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newy[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newz[domain[[i,1]],domain[[i,2]],domain[[i,3]]]},{i,1,Length[domain]}];
Clear[partimage]
partimage=Table[image[[i]],{i,1,10}];
PrependTo[partimage,{"x coordinates","y coordinates","z coordinates"}];
Print["The first ten image points are: "]
Print[TableForm[partimage]]

Let's suppose that there is a tumor in the region bounded approximately by z = 100 -0.5(Sqrt[(x^2)+(y^2)]) and z = .005((x^2)+(y^2)).  We will first plot the object in three-space and then separate out the image points that have been mapped to that region and plot them in green.

Clear[r]
p1=RevolutionPlot3D[100-0.5` r,{r,0,100},{θ,0,2 π},DisplayFunction->Identity];
p2=RevolutionPlot3D[0.005` (r^2),{r,0,100},{θ,0,2 π},DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction,ViewPoint->{-0.637`,-0.917`,-0.018`}]
listout={};
listin={};
Do[If[domain[[i,3]]<100-0.5` √(((domain[[i,1]])^2)+((domain[[i,2]])^2))&&domain[[i,3]]>0.005` (((domain[[i,1]])^2)+((domain[[i,2]])^2)),AppendTo[listin,{image[[i,1]],image[[i,2]],image[[i,3]]}],AppendTo[listout,{image[[i,1]],image[[i,2]],image[[i,3]]}]],{i,1,Length[domain]}]
pout=ListPointPlot3D[listout,DisplayFunction->Identity];
pin=ListPointPlot3D[listin,PlotStyle->{RGBColor[0,1,0],PointSize[0.02`]},DisplayFunction->Identity];
Show[pin,pout,DisplayFunction->$DisplayFunction]

Suppose, for example, that a spot is found at the image point T1 (10,20,30). What does that tell you about the location of the possible growth? Based on this image, you cannot precisely identify the placement of the tumor, because this mapping is NOT one-to-one.

maptoT1=Solve[{newx[x,y,z]==10,newy[x,y,z]==20,newz[x,y,z]==30},{x,y,z}]

Since this result gives you some information, but not enough, you must take another X-ray from a different angle. 
We will set up a mapping from a different perspective. In the example given below, we plot the new perspective and verify that the points all map to a plane.

Clear[x,y,z,newx,newy,newz,domain,image,listin,listout,pin,pout]
newx[x_,y_,z_]:=-x+y
newy[x_,y_,z_]:=2 x-3 y+2 z
newz[x_,y_,z_]:=-y+2 z
Print["The number of domain and image points is ",numberofpoints=1000]
domain=Table[RandomInteger[{-100,100},3],{i,1,numberofpoints}];
image=Table[{newx[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newy[domain[[i,1]],domain[[i,2]],domain[[i,3]]],newz[domain[[i,1]],domain[[i,2]],domain[[i,3]]]},{i,1,Length[domain]}];
listout={};
listin={};
Do[If[domain[[i,3]]<100-0.5` √(((domain[[i,1]])^2)+((domain[[i,2]])^2))&&domain[[i,3]]>0.005` (((domain[[i,1]])^2)+((domain[[i,2]])^2)),AppendTo[listin,{image[[i,1]],image[[i,2]],image[[i,3]]}],AppendTo[listout,{image[[i,1]],image[[i,2]],image[[i,3]]}]],{i,1,Length[domain]}]
pout=ListPointPlot3D[listout,DisplayFunction->Identity];
pin=ListPointPlot3D[listin,PlotStyle->{RGBColor[0,1,0],PointSize[0.02`]},DisplayFunction->Identity];
Show[pin,pout,DisplayFunction->$DisplayFunction]

If you were somehow able to detect that the SAME spot now shows up on the image at T2(-90,250,70), would that tell you about the possible location of the tumor? As before, we can get information on the point in the domain that mapped to this image point.

maptoT2=Solve[{newx[x,y,z]==-90,newy[x,y,z]==250,newz[x,y,z]==70},{x,y,z}]

If we put this information together with the previous information, we can locate the point on the organ identified. Essentially, we are finding the intersection of two straight lines in three-dimensional space.

var=Solve[{maptoT1[[1,1,2]]==maptoT2[[1,1,2]],maptoT1[[1,2,2]]==maptoT2[[1,2,2]]},z]

Let's put this result for z into the specifications for the first line we found.

Print["x = ",(maptoT1[[1,1,2]]/.var)[[1]]]
Print["y = ",(maptoT1[[1,2,2]]/.var)[[1]]]
Print["z = ",var[[1,1,2]]]

This specifies one of the points in the three-dimensional object where the tumor is present.

You should notice in the analysis above how difficult it is to identify the point in the second mapping that is the image of the SAME point as the one we first found. Because of this problem, many more than simply two X-rays would need to be taken to determine the precise location of a tumor The technology that makes this process feasible is called a CT-Scan, commonly referred to as a CAT-Scan. The procedure above reflects the fundamental mathematics used to interpret CAT-Scans.

Part V: Linear Algebra Approach (optional; not in calculus text)

\[MathematicaIcon] About Mathematica

Clear[m,eig] 
m={{1,1,2},{2,1,-1},{3,2,1}};MatrixForm[m] 
eig=Eigensystem[m];Print["The eigenvalues and corresponding eigenvectors for the given matrix are: ",N[eig]]
image={}; Do[{new=m.RandomInteger[{-100,100},3],AppendTo[image,new]},{1000}]

Note the presence of the 0 eigenvalue.

The kernel of a linear transformation is the set of elements in the domain that map to the 0 element (0, 0, 0) in the image.

Clear[a,b,c]
kernel={a,b,c};
solk=Solve[m.kernel=={0,0,0},{a,b,c}]

If c = 1, note how this vector compares to the eigenvector associated with the eigenvalue of 0.

c=1;
kernel/.solk

You can find the normal to the plane formed by the image by taking the cross product of the eigenvectors associated with the nonzero eigenvalues. In this case, these are the first and second eigenvectors.

norm=Cross[eig[[2,1]],eig[[2,2]]];
scalednorm=norm/norm[[1]]//Simplify

Check to see if all the points in the image lie in the plane.

works=0;
Do[If[scalednorm[[1]]image[[i,1]]+scalednorm[[2]]image[[i,2]]+scalednorm[[3]]image[[i,3]]==0,works=works+1],{i,1,1000}] 
works

What do you suppose the existence of a 0 eigenvalue has to do with the fact that the image points all lie in a plane, even though the domain consists of points in three dimensions?

\[MathematicaIcon] About Mathematica

If you have studied linear algebra, you can learn in this section some of the behind-the-scenes mathematics in the study of linear transformations and how Mathematica can help you with this analysis.
The command Eigensystem[m] in Mathematica gives a list consisting of two parts. The first part contains the eigenvalues of the matrix for m and the second part gives the set of eigenvectors for the eigenvalues in corresponding order. 
 Go back.