Mathematica Code for: “Transformations and Projections in Computer Graphics,” 2006, by David  Salomon.

Code Examples:

(* Sine Values as 16-Bit Integers. Page 23. *)
Clear[d2r]
d2r=Pi/180;
s1=Table[Round[N[16384*Sin[i*d2r]]],{i,0,30}];
a1=Table[i,{i,0,30}];
s2=Table[Round[N[16384*Sin[i*d2r]]],{i,31,60}];
a2=Table[i,{i,31,60}];
s3=Table[Round[N[16384*Sin[i*d2r]]],{i,61,90}];
a3=Table[i,{i,61,90}];
tab1=MatrixForm[Transpose[{Prepend[a1,"θ"],Prepend[s1,"Sin θ"]}]];
tab2=MatrixForm[Transpose[{Prepend[a2,"θ"],Prepend[s2,"Sin θ"]}]];
tab3=MatrixForm[Transpose[{Prepend[a3,"θ"],Prepend[s3,"Sin θ"]}]];
Quiet[StringJoin[{tab1},{tab2},{tab3}]];
Framed[%,Background->LightCyan];
Labeled[%,"Sine Values as 16-Bit Integers",{Top}];
Framed[%,Background->LightYellow]

(* CORDIC Rotation. Fig. 1.12, page 26. *)
N[Table[Product[(2^(-2i)+1)^(-1/2),{i,0,n}],{n,0,16}],17];
Column[%,Alignment->"."];
Framed[%,Background->LightBlue];
Labeled[%,"CORDIC Rotation",{Top}];
Framed[%,Background->LightOrange,
FrameStyle->Directive[Thick,Brown]]

(* CORDIC Rotations. Fig.1.13, page 28. *)
Clear[t,d,x,y,z]
t=Table[ArcTan[2.^{-i}],{i,0,15}];(*Arctans in radians.*)
d=1;x=2.1;y=0.34;z=46. Degree;
Do[{Print[i,",",x,",",y,",",z,",",d],
xn=x+y d 2^{-i},yn=y-x d 2^{-i},
zn=z-d t[[i+1]],d=Sign[zn],x=xn,y=yn,z=zn},{i,0,14}]
Print[0.60725x,",",0.60725y]

(*Mathematica Code for a General Rotation. Fig. 127, page 50.*)
Clear[tm, a, b, c, t, x, y, z]
tm = Sqrt[x^2 + y^2];
a = {{x/tm, -y/tm, 0}, {y/tm, x/tm, 0}, {0, 0, 1}};
b = {{z, 0, Sqrt[1 - z^2]}, {0, 1, 0}, {-Sqrt[1 - z^2], 0, z}};
c = {{Cos[t], -Sin[t], 0}, {Sin[t], Cos[t], 0}, {0, 0, 1}};
FullSimplify[a.b.c.Transpose[b].Transpose[a] /. x^2 + y^2 -> 1 - z^2]

(*Computing Three Givens Matrices.Fig.1.28,page 51.(Mathematica code.)*)
n=3;
a={{.5774,-.5774,-.5774},{.5774,.7886,-.2115},{.5774,-.2115,.7886}};
(* Rotation from 1,1,1 to x-axis *)
qq=IdentityMatrix[n];
Do[
Do[
Tt=IdentityMatrix[n];
      dd=Sqrt[a[[j,j]]^2+a[[i,j]]^2];
      coss=a[[j,j]]/dd;ssin=a[[i,j]]/dd;
Tt[[j,j]]=coss;Tt[[j,i]]=ssin;
Tt[[i,j]]=-ssin;Tt[[i,i]]=coss;
a=Tt.a;
       qq=qq.Transpose[Tt],
  { i,n,j+1,-1}],
{ j,1,n-1}];
qq
a

Chop[%]

(*Rotating point (1,1,1) to the x axis.Fig.1.29,page 51.(Mathematica code.)*)
T1={{0.7071,0,0.7071},{0,1,0},{-0.7071,0,0.7071}};
T2={{0.8165,0.5774,0},{-0.5774,0.8165,0},{0,0,1}};
T3={{1,0,0},{0,0.9660,0.2587},{0,-0.2587,0.9660}};
p={1,1,1};
a=T1.p
b=T2.a
c=T3.b

Two examples using Matlab code:


Computing Three Givens Matrices. Fig. 1.28, page 51. (Matlab code)

n=3;
A=[.5774,-.5774,-.5774; .5774,.7886,-.2115; .5774,-.2115,.7886]
% Rotation from 1,1,1 to x-axis
Q=eye(n);
for j=1:n-1,
for i=n:-1:j+1,
T=eye(n);
D=sqrt(A(j,j)^2+A(i,j)^2);
cos=A(j,j)/D; sin=A(i,j)/D;
T(j,j)=cos; T(j,i)=sin; T(i,j)=-sin; T(i,i)=cos; T
A=T*A;
Q=Q*T';
end;
end;
Q
A


Rotating point (1,1,1) to the x axis. Fig. 1.29, page 51. (Matlab code)

T1=[0.7071,0,0.7071; 0,1,0; -0.7071,0,0.7071];
T2=[0.8165,0.5774,0; -0.5774,0.8165,0; 0,0,1];
T3=[1,0,0; 0,0.9660,0.2587; 0,-0.2587,0.9660];
p=[1;1;1];
a=T1*p
b=T2*a
c=T3*b


(Note: The rest of this file is done using Mathematica code.)

{0,1,-4,1}.(({{1, 0, 0, 0}, {0, 0, 1, 0}, {0, -1, 0, 0}, {0, 0, 0, 1}})).(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 3, 1}})).(({{1, 0, 0, 0}, {0, 1/2, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}})).(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 1}}))(* Page 94, Mathematica code. *)

(* tests of matrix T_g, equation {eq:pers.gen5}. Code begins on page 110. *)
Clear[k,r,a,b,c,d,e,f,T]
k=3.; r=1/k;
{a,b,c}={0,0,-k}; {d,e,f}=Normalize[{-1,1,k}]
T={{(e^2+f+f^2)/(1+f),-d e/(1+f),0,d r}, {-d e/(1+f),(d^2+f+f^2)/(1+f),0,e r},{-d,-e,0,f r},
{(c d+b d e-a e^2-a f+c d f-a f^2)/(1+f),(-b d^2+c e+a d e-b f+c e f-b f^2)/(1+f),0,-(a d+b e+c f)  r}};
{-1,1,0,1}.T
{a,b,c}={0,2k,-2k}; {d,e,f}=Normalize[{2k,-2k,2k}]
T={{(e^2+f+f^2)/(1+f),-d e/(1+f),0,d r}, {-d e/(1+f),(d^2+f+f^2)/(1+f),0,e r},{-d,-e,0,f r},
{(c d+b d e-a e^2-a f+c d f-a f^2)/(1+f),(-b d^2+c e+a d e-b f+c e f-b f^2)/(1+f),0,-(a d+b e+c f)  r}};
{2k,0,0,1}.T
(* for exercise that follows the tests of matrix T_g *)
{a,b,c}={0,2k,-k}; {d,e,f}=Normalize[{0,-1,-1}]
T={{(e^2+f+f^2)/(1+f),-d e/(1+f),0,d r}, {-d e/(1+f),(d^2+f+f^2)/(1+f),0,e r},{-d,-e,0,f r},
{(c d+b d e-a e^2-a f+c d f-a f^2)/(1+f),(-b d^2+c e+a d e-b f+c e f-b f^2)/(1+f),0,-(a d+b e+c f)  r}};
{0,0,-4k,1}.T

Chop[Out[51]]

Chop[Out[54]]

N[Out[55]]

(* Used in the discussion of the top vector in general perspective. (Ref: 3.5 General Perspective, page 96.) Fig.3.35, pages 112-113*)
Clear[a,h,T,i]

a=1/(Sqrt[2]);
h=a; i=a;T=(({{i, h, 0, 0}, {-h, i, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}));
  {h,i,0,1}.T

(* Code to check matrix T_g for the case 1+f=0. Page 112. *)
Clear[a,b,c,r,k,x,y,z,T];
r=1/k; 
{a,b,c}={1,1,1};
T={{-1,0,0,0},{0,-1,0,0},{0,0,0,-r},{a,b,0,c r}};
{1,1,-1,1}.T

(* For exercise that follows the tests of matrix T_g *)
Clear[k,r,a,b,c,d,e,f,T]
{a,b,c}={0,2k,-k};
{d,e,f}=Normalize[{0,-1,-1}]
T={{(e^2+f+f^2)/(1+f),-d e/(1+f),0,d r}, {-d e/(1+f),(d^2+f+f^2)/(1+f),0,e r},{-d,-e,0,f r},
{(c d+b d e-a e^2-a f+c d f-a f^2)/(1+f),(-b d^2+c e+a d e-b f+c e f-b f^2)/(1+f),0,-(a d+b e+c f)  r}};
{0,0,-4k,1}.T

(* Code for page 113. *)
Clear[a,h,i,T]
a=1/Sqrt[2];h=a;i=a;
T={{i,h,0,0},{-h,i,0,0},{0,0,1,0},{0,0,0,1}};
{h,i,0,1}.T

Drop[%,-1]

Append[%,1]

(* Display two cubes as a stereo pair. Fig.3.42(a), page 126. *)
Clear[Trg,Tlf,pt,e,r,qt];
(*Tlf={{1,0,0,0},{0,1,0,0},{0,0,0,r},{e,0,0,1}};*)
Trg={{1,0,0,0},{0,1,0,0},{0,0,0,r},{-e,0,0,1}}; 
pt={{1,1,1,1},{-1,1,1,1},{1,-1,1,1},{-1,-1,1,1},{1,1,-1,1},{-1,1,-1,1},{1,-1,-1,1},{-1,-1,-1,1},{1,1,1,1}};
e=.1;r=3;
qt=Table[0,{i,9},{j,4}];
Do[qt[[i]]=pt[[i]].Trg,{i,1,9}]; (* Use Tlf for other image. *)
Do[qt[[i,1]]=qt[[i,1]]/qt[[i,4]],{i,1,9}];
Do[qt[[i,2]]=qt[[i,2]]/qt[[i,4]],{i,1,9}];
ListPlot[Table[{qt[[i,1]],qt[[i,2]]},{i,1,9}],Joined->True,AspectRatio->1,ColorFunction->"ArmyColors",PlotStyle->Thick,Background->Lighter[LightGreen],Frame->True,Axes->False,FrameLabel->"Fig. 3.42(a)",FrameTicks->False]

(* Display two cubes as a stereo pair. Fig.3.42(b), page 126. *)
Clear[Trg,Tlf,pt,e,r,qt];
Tlf={{1,0,0,0},{0,1,0,0},{0,0,0,r},{e,0,0,1}};
(*Trg={{1,0,0,0},{0,1,0,0},{0,0,0,r},{-e,0,0,1}}; *)
pt={{1,1,1,1},{-1,1,1,1},{1,-1,1,1},{-1,-1,1,1},{1,1,-1,1},{-1,1,-1,1},{1,-1,-1,1},{-1,-1,-1,1},{1,1,1,1}};
e=.1;r=3;
qt=Table[0,{i,9},{j,4}];
Do[qt[[i]]=pt[[i]].Tlf,{i,1,9}]; (* Use Trg for other image *)
Do[qt[[i,1]]=qt[[i,1]]/qt[[i,4]],{i,1,9}];
Do[qt[[i,2]]=qt[[i,2]]/qt[[i,4]],{i,1,9}];
ListPlot[Table[{qt[[i,1]],qt[[i,2]]},{i,1,9}],Joined->True,AspectRatio->1,ColorFunction->"DarkRainbow",PlotStyle->Thick,Background->Lighter[LightCyan],Frame->True,Axes->False,FrameLabel->"Fig. 3.42(b)",FrameTicks->False]

(* False perspective, nonlinear projections. Pages 145-146, Fig.4.1. *)
Clear[a,b,tabl1,tabl2,f];
a=.001; b=.001;
tabl1=Table[Point[{i,j}], {i,-10,10, 5}, {j,-10,10, 5}];
f[x_,y_]:=-(1/2.)Exp[-a x^2-b y^2];
tabl2=Table[Point[{i,j}/(1+f[i,j])], {i,-10,10, 5}, {j,-10,10, 5}];
Show[Graphics[{PointSize[0.02],Black,tabl1}],Background->Lighter[LightBlue],GridLines->Automatic,Frame->True,FrameTicks->False]
Show[Graphics[{PointSize[0.02],Red,tabl2}],Background->LightYellow,GridLines->Automatic,Frame->True,FrameLabel->"Fig. 4.1: Moving Points in False Perspective.",FrameTicks->False]

(* Hemispherical Fisheye Projection. Page 149, Fig. 4.3. *)
Clear[k,n,P,Q,L]
k=10; n=50;
scal[q_]:=(k Tan[ArcTan[q/k]/2])/q;
P=Table[{Random[Real,{-10.,10.}],Random[Real,{-10.,10.}]}, {n}];
Q=Table[Sqrt[P[[i]].P[[i]]], {i,n}];
L=Table[{Hue[i/20],Line[{P[[i]],scal[Q[[i]]] P[[i]]}]}, {i,n}];
Show[Graphics[L], Graphics[Circle[{0,0},10]],Graphics[Point[{0,0}]], AspectRatio->1,Background->Lighter[Lighter[LightBlue]]]
Export["xyz.pdf",%](*This line makes a pdf file named "xyz.pdf" for the graph below.*)

(* Animation: Hemispherical Fisheye Projection. Page 149, Fig. 4.3. *)
Animate[Clear[k,n,P,Q,L]
k=10; n=50;
scal[q_]:=(k Tan[ArcTan[q/k]/2])/q;
P=Table[{Random[Real,{-10.,10.}],Random[Real,{-10.,10.}]}, {n}];
Q=Table[Sqrt[P[[i]].P[[i]]], {i,n}];
L=Table[{Hue[i/20],Line[{P[[i]],scal[Q[[i]]] P[[i]]}]}, {i,n}];
tx=Graphics[{Text[Style["Hemispherical Fisheye Projection",FontFamily->"Jokewood",FontSize->14,Blue,Bold],{0,-11}]}];
Show[Graphics[tx],Graphics[L], Graphics[Circle[{0,0},10]],Graphics[Point[{0,0}]], AspectRatio->1,Background->Lighter[Lighter[LightBlue]]],{k,0,6}]

(*Code for Table 4.6. Page 151*)
Clear[k,c1,c2,c3,c4,angl,data]
k = 10;
c1=Table[i,{i,0,4}];
c2=angl = {0,22.5, 45., 67.5, 89.};
c3=k Tan[angl Degree];
c4=k Tan[angl/2 Degree];
data=Style[MatrixForm[Transpose[{Prepend[c1,"#"],Prepend[c2,"θ"],Prepend[c3,"(^R)n"],Prepend[c4,"(^r)n"]}]],{16,Darker[Purple]}];
Framed[data,Background->Lighter[LightPurple]];
Labeled[%,"Table 4.6",{Top}];
Framed[%,Background->Lighter[Lighter[LightGray]]]

(* Offaxis Hemispherical Fisheye Projection. Page 157, Fig. 4.12(b). *)
Clear[k,n,P,Q,L,Pt];
k=10; n=50;
scal[q_]:=(k Tan[ArcTan[q/k]/2])/q;
P=Table[{Random[Real,{-10.,10.}],Random[Real,{-10.,10.}]}, {n}];
x=-5; y=5; (* Location of viewer *)
Pt=P-Table[{x,y}, {n}];
Q=Table[Sqrt[Pt[[i]].Pt[[i]]], {i,n}];
L=Table[{Red,Line[{P[[i]]+{x,y},(scal[Q[[i]]] P[[i]])+{x,y}}]}, {i,n}];
Show[Graphics[L], Graphics[Circle[{0,0},k]],
Graphics[{AbsolutePointSize[5],Point[{0,0}]}],
Graphics[{AbsolutePointSize[5],Point[{x,y}]}],
 AspectRatio->Automatic, PlotRange->All,Background->LightPurple]

4.3.  Circle Inversion

Go to the Wolfram Demonstrations Project and see Orthogonal Circle Inversion. This is a free user interactive downloadable Mathematica notebook.
Just click on the hyperlink below.

Hyperlink["http://demonstrations.wolfram.com/OrthogonalCircleInversion/"]

(* Telescopic Projection. Fig. 4.41(a) Twenty Nonuniformly Spaced Points. Page 195. *) 
Clear[k,z,l,r];
k=10.; 
Table[k z/(z+k), {z,0,100,5}]
Table[%[[i+1]]-%[[i]], {i,1,20}]
Table[{Purple,Point[{%%[[i]],0}]}, {i,1,21}];
Show[Graphics[%],Background->LightYellow,PlotLabel->"Twenty Nonuniformly Spaced Points",Frame->True]

l=20.;r=0.1;(*Varing Heights. Fig. 4.41(b), Page 196.*)
Table[l(1-(z r/(z+l))), {z,0,100,5}]
Table[%[[i]]-%[[i+1]], {i,1,20}]
Table[Line[{{i,17},{i,%%[[i]]}}], {i,1,21}]
Show[Graphics[%],Background->LightCyan,PlotLabel->"Varing Heights",Frame->True]

(*1.27(a), page 251. 60 Degrees.*)
t14=2^14;
Print["(x*=",(8192-(2*14189.))/t14,",y*=",(14189.+(2*8192))/t14,")"]
Print["(x*=",Cos[60Degree]-2.Sin[60Degree],
",y*=",Sin[60Degree]+2.Cos[60Degree],")"]

(*1.27(b), page 251-252. 80 Degrees.*)
t14=2^14;
Print["(x*=",(2845-(2*16135.))/t14,",y*=",(16135.+(2*2845))/t14,")"]
Print["(x*=",Cos[80Degree]-2.Sin[80Degree],
",y*=",Sin[80Degree]+2.Cos[80Degree],")"]

(* Code for exercise {ex:pers.sqrt8}. Page 264 *)
Clear[a,b,c,d,e,f,r,T]
{a,b,c}={0,1.,0}; 
{d,e,f}=Normalize[{0,1,1}]
T={{(e^2+f+f^2)/(1+f),-d e/(1+f),0,d r}, {-d e/(1+f),(d^2+f+f^2)/(1+f),0,e r},{-d,-e,0,f r},
{(c d+b d e-a e^2-a f+c d f-a f^2)/(1+f),(-b d^2+c e+a d e-b f+c e f-b f^2)/(1+f),0,-(a d+b e+c f)  r}};
{0,1,10,1}.T

(* Exercise to compute matrix \bfT_{123}. Answers 3.25, pages 264-265. *)
Clear[a,b,c,d,e,f,k,r,k]
k=Sqrt[2]; r=1/k;
{a,b,c}={0,2k,-2k}; {d,e,f}={0,-1,-.99};
T123=
(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {-a, -b, -c, 1}})).(({{(e^2+f-f^3-e^2 f)/(1-f^2), -e d/(1+f), d, 0}, {-e d/(1+f), (d^2+f-f^3-d^2 f)/(1-f^2), e, 0}, {-d, -e, f, 0}, {0, 0, 0, 1}})).(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, r}, {0, 0, 0, 1}}))
{0,2k-3,-2k+1,1}.T123

(* Exercise for Hemispherical Fisheye Projection. (4.1. Pages 267-268). *)
Clear[k,q,scal]
k=1;
scal[q_]:=(k Tan[ArcTan[q/k]/2])/q;
{scal[1.],scal[10.], scal[100.], scal[1000.], scal[10000.]}

For further information about books by David Salomon, and free downloadable Mathematica notebooks and other information,
just click on the hyperlink below.

Hyperlink["http://www.davidsalomon.name/"]

This notebook was done with Mathematica 6.0.1.  10-28-2008. 10-29-2008 and 10-30-2008.