3Dグラフィックスではレンダリングパイプラインというものがあります。
最近のopenglなどは固定のパイプラインを使わずに
プログラマが自分でシェーダを記述して
コーディングする形が主流になっていますが
ここでは3Dグラフィックの理解のために
maximaで3Dのレンダリングパイプラインを計算するコードを
書いてみましょう。

3Dのレンダリングパイプライン
/* 補助関数 */
fold(f,init,lst):=if lst=[] then init else fold(f,f(init,first(lst)),rest(lst));
take(n,l):=if n<=0 then [] else cons(first(l),take(n-1,rest(l)));
/* 正規化 */
norm(x):=sqrt(fold(lambda([a,x],a+x*x),0,x));
/* クロス積 */
cross( u, v ):= [u[2]*v[3] - u[3]*v[2], u[3]*v[1] - u[1]*v[3], u[1]*v[2] - u[2]*v[1] ];

/* 斉次座標の回転や移動、スケーリング*/
rotatex(x):=matrix([1,0,0,0],[0,cos(x),sin(x),0],[0,-sin(x),cos(x),0],[0,0,0,1]);
rotatey(x):=matrix([cos(x),0,-sin(x),0],[0,1,0,0],[sin(x),0,cos(x),0],[0,0,0,1]);
rotatez(x):=matrix([cos(x),sin(x),0,0],[-sin(x),cos(x),0,0],[0,0,1,0],[0,0,0,1]);
move(x,y,z):=matrix([1,0,0,0],[0,1,0,0],[0,0,1,0],[x,y,z,1]);
scale(x,y,z):=matrix([x,0,0,0],[0,y,0,0],[0,0,z,0],[0,0,0,1]);

rotate_matrix(c,t):=block([x,y,z],z:t-c,z:z/norm(z),x:cross(z,[0,0,1]),x:x/norm(x),y:cross(z,x),x:endcons(-(x.c),x),y:endcons(-(y.c),y),z:endcons(-(z.c),z),transpose(matrix(x,y,z,[0,0,0,1])));
perspective_projection(x):=matrix([1,0,0,0],[0,1,0,0],[0,0,0,-1/norm(x)],[0,0,0,1]);
viewport(w,h):=scale(w/2,-h/2,1) . move(w/2,h/2,0);
mat3do(m,c,at,width,height):=block([model,camera,mat,w],model:matrix(endcons(1,m)),camera:matrix(c),mat:model . rotate_matrix(c,at) . perspective_projection(camera[1]) . viewport(width,height),w:mat[1][4],mat/w);
mat3dob(m,c,at,w,h):=block([r],r:mat3do(m,c,at,w,h),[float(r[1][1]),float(r[1][2])]);

fpprintprec:5;
mat3dob([1,1,1],[0,0,-1],[2,2,1],100,100);



画像処理の教科書では、計算済みの複雑な行列のみが載っていたりしますが
自分でコードを書いてみると、行列を何度か掛け合わせただけということが
わかります。