在計算機圖形學中，矩陣乘法有着很好的應用。圖形的變換可以通過構造相應的矩陣進行計算來完成。
我們知道，平面上的元素，就是點、線、面，而線就是由一個個點組成的，面是由一條條線組成的，所以歸根結底，平面上所有的圖形都是由點組成的。在坐標系中，一個點就是由一對x，y值組成的，p = {x, y}。在平面上，過兩點間的，可以畫一條直線，所以我們一般通過 兩個點p1, p2可定義一條直線，e = {p1, p2}，而圖形呢，則是由眾多的點和點之間的的線段組成的。所以，平面上的圖形變換，本質是點坐標位置的變換。
在平面中，常用的基本圖形變換有以下四種：
1）平移（Translation）。
設點(x,y)水平向右平移dx個單位，垂直向上平移dy個單位。表示點（x,y）的矩陣P和構造的平移變換矩陣T如下：
則點(x,y)的平移后的新位置P'可以通過矩陣乘法計算出來。
2）縮放（Scale）。
設點(x,y)在水平方向和垂直方向的縮放比例分別為Sx和Sy。表示點（x,y）的矩陣P和構造的縮放變換矩陣S如下：
則點(x,y)的縮放后的新位置P'可以通過矩陣乘法計算出來。
3）旋轉（Rotation）。
設點(x,y)繞原點逆時針旋轉α角度。表示點（x,y）的矩陣P和構造的旋轉變換矩陣R如下：
旋轉矩陣的構造原理參見下圖。
則點(x,y)繞原點逆時針旋轉后的新位置P'可以通過矩陣乘法計算出來。
4）反射（Reflect）。
反射變換分為關於原點反射、關於x軸反射和關於y軸反射三種。
點 p(x,y) 關於原點反射后得到的點為po'(-x,-y)。
點 p(x,y) 關於x軸反射后得到的點為px'(x,-y)，也稱為上下翻轉。
點 p(x,y) 關於x軸反射后得到的點為py'(-x,y)，也稱為左右翻轉。
三種反射的矩陣計算如下：
【例1】點的變換。
描述
平面上有不超過10000個點，坐標都是已知的，現在可能對所有的點做以下幾種操作：
平移一定距離(M)，相對X軸上下翻轉(X)，相對Y軸左右翻轉(Y)，坐標縮小或放大一定的倍數(S)，所有點對坐標原點逆時針旋轉一定角度(R)。
操作的次數不超過10000次，求最終所有點的坐標。
輸入
測試數據的第一行是兩個整數N,M，分別表示點的個數與操作的個數(N,M<=10000)
隨后的一行有N對數對，每個數對的第一個數表示一個點的x坐標，第二個數表示y坐標，這些點初始坐標大小絕對值不超過100。
隨后的M行，每行代表一種操作，行首是一個字符：
首字符如果是M,則表示平移操作，該行后面將跟兩個數x,y，表示把所有點按向量(x,y)平移;
首字符如果是X，則表示把所有點相對於X軸進行上下翻轉;
首字符如果是Y，則表示把所有點相對於Y軸進行左右翻轉;
首字符如果是S，則隨后將跟一個數P,表示坐標放大P倍;
首字符如果是R，則隨后將跟一個數A,表示所有點相對坐標原點逆時針旋轉一定的角度A(單位是度)
輸出
每行輸出兩個數，表示一個點的坐標(對結果四舍五入到小數點后1位，輸出一位小數位）
點的輸出順序應與輸入順序保持一致
樣例輸入
2 5
1.0 2.0 2.0 3.0
X
Y
M 2.0 3.0
S 2.0
R 180
樣例輸出
-2.0 -2.0
0.0 0.0
（1）編程思路。
直接按前面的介紹內容，分別構造5類變換矩陣。定義矩陣ans初始為單位矩陣。然后每輸入m種操作中的一種時，根據操作類型構造相應的變換矩陣temp，然后ans右乘temp，即ans=temp*ans。輸入完m種操作后，得到的ans矩陣就是每個點的變換矩陣，進行變換矩陣與點的矩陣乘法就可以得到每個點的最終坐標。
（2）源程序。
#include <stdio.h>
#include <string.h>
#include <math.h>
#define PI acos(-1.0)
#define MAXN 10005
struct Matrix
{
double mat[4][4]; // 存儲矩陣中各元素
};
struct Point
{
double x ,y ;
};
Matrix matMul(Matrix a ,Matrix b,int n)
{
Matrix c;
memset(c.mat,0,sizeof(c.mat));
int i,j,k;
for (k = 1; k<=n ; k++)
for (i=1 ;i<=n ; i++)
if (a.mat[i][k]!=0)
for (j = 1 ;j<=n ;j++)
c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) ;
return c;
}
int main()
{
char ch ;
Point p[MAXN];
Matrix ans,temp;
double a ;
int n,m,i,j;
scanf("%d%d",&n,&m);
for (i=0 ;i<n;i++)
scanf("%lf%lf" ,&p[i].x ,&p[i].y);
memset(ans.mat,0,sizeof(ans.mat)); // 初始為單位陣
ans.mat[1][1]=ans.mat[2][2]=ans.mat[3][3]=1;
for (i=0 ;i<m;i++)
{
memset(temp.mat,0,sizeof(temp.mat)); // 初始為單位陣
temp.mat[1][1]=temp.mat[2][2]=temp.mat[3][3]=1;
getchar();
scanf("%c",&ch);
if (ch == 'M') // 平移
scanf("%lf%lf" ,&temp.mat[1][3] ,&temp.mat[2][3]);
else if (ch == 'X') // 相對X軸上下翻轉
temp.mat[2][2]=-1;
else if (ch == 'Y') // 相對於Y軸左右翻轉
temp.mat[1][1] = -1 ;
else if (ch == 'S') // 坐標放大P倍
{
scanf("%lf" ,&temp.mat[1][1]) ;
temp.mat[2][2] = temp.mat[1][1] ;
}
else if (ch == 'R') // 相對坐標原點逆時針旋轉一定的角度A
{
scanf("%lf" ,&a) ;
a = (a*PI)/180.0 ;
temp.mat[1][1] = temp.mat[2][2] = cos(a) ;
temp.mat[1][2] = -sin(a) ;
temp.mat[2][1] = sin(a) ;
}
ans=matMul(temp,ans,3);
}
for (i=0 ;i<n ; i++)
{
for (j=1 ;j<=2; j++)
{
if (j!=1) printf(" ");
printf("%.1lf" ,(ans.mat[j][1]*p[i].x+ans.mat[j][2]*p[i].y+ans.mat[j][3]));
}
printf("\n");
}
return 0 ;
}
【例2】openGL （HDU 3320）。
Problem Description
Jiaoshou selected a course about “openGL” this semester. He was quite interested in modelview, which is a part of “openGL”. Just using three functions, it could make the model to move, rotate and largen or lessen. But he was puzzled with the theory of the modelview. He didn’t know a vertex after several transformations where it will be.
Now, He tells you the position of the vertex and the transformations. Please help Jiaoshou find the position of the vertex after several transformations.
Input
The input will start with a line giving the number of test cases, T.
Each case will always begin with “glBegin(GL_POINTS);”.Then the case will be followed by 5 kinds of function.
1. glTranslatef(x,y,z);
This function will translate the vertex(x’,y’,z’) to vertex(x+x’,y+y’,z+z’).
2. glRotatef(angle,x,y,z);
This function will turn angle radians counterclockwise around the axis (0,0,0)->(x,y,z).
3. glScalef(x,y,z);
This function wiil translate the vertex(x’,y’,z’) to vertex(x*x’,y*y’,z*z’).
4. glVertex3f(x,y,z);
This function will draw an initial vertex at the position(x,y,z). It will only appear once in one case just before “glEnd();”. In openGL, the transformation matrices are right multiplied by vertex matrix. So you should do the transformations in the reverse order.
5. glEnd();
This function tells you the end of the case.
In this problem angle,x,y,z are real numbers and range from -50.0 to 50.0. And the number of functions in each case will not exceed 100.
Output
For each case, please output the position of the vertex after several transformations x,y,z in one line rounded to 1 digits after the decimal point , separated with a single space. We guarantee that x,y,z are not very large.
Sample Input
1
glBegin(GL_POINTS);
glScalef(2.0,0.5,3.0);
glTranslatef(0.0,1.0,0.0);
glVertex3f(1.0,1.0,1.0);
glEnd();
Sample Output
2.0 1.0 3.0
Hint
In this sample, we first let the vertex do “glTranslatef(x,y,z);” this function, then do “glScalef(x,y,z)”.
（1）編程思路。
題目的意思是：給出三維空間一個點的坐標，對它進行多種變換（平移、縮放、繞過原點的任意軸旋轉），輸出它的最終坐標。需要注意的是它給出各變換的操作順序是反過來的（見題目Hint）。
根據輸入構造好三維坐標點的平移、縮放和旋轉轉換矩陣，用矩陣乘法乘起來就是最終坐標。
（2）源程序。
#include <stdio.h>
#include <math.h>
#include <string.h>
int cnt;
char opstr[105][105];
struct Matrix
{
double mat[5][5]; // 存儲矩陣中各元素
};
Matrix matMul(Matrix a ,Matrix b,int n)
{
Matrix c;
memset(c.mat,0,sizeof(c.mat));
int i,j,k;
for (k = 1; k<=n ; k++)
for (i=1 ;i<=n ; i++)
if (a.mat[i][k]!=0)
for (j = 1 ;j<=n ;j++)
c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]);
return c;
}
Matrix translate(Matrix p, double x, double y, double z) // 平移
{
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
for (int i=1;i<=4;i++)
tmp.mat[i][i]=1;
tmp.mat[1][4] = x;
tmp.mat[2][4] = y;
tmp.mat[3][4] = z;
p=matMul(tmp,p,4);
return p;
}
Matrix scale(Matrix p, double a, double b, double c) //縮放
{
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
for (int i=1;i<=4;i++)
tmp.mat[i][i]=1;
tmp.mat[1][1] = a;
tmp.mat[2][2] = b;
tmp.mat[3][3] = c;
p=matMul(tmp,p,4);
return p;
}
Matrix rotate(Matrix p, double x, double y, double z, double d) // 旋轉
{
double di = sqrt(x*x+y*y+z*z); // 單位化
x/=di, y/=di, z/=di;
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
tmp.mat[1][1] = (1-cos(d))*x*x + cos(d);
tmp.mat[1][2] = (1-cos(d))*x*y - z*sin(d);
tmp.mat[1][3] = (1-cos(d))*x*z + y*sin(d);
tmp.mat[2][1] = (1-cos(d))*x*y + z*sin(d);
tmp.mat[2][2] = (1-cos(d))*y*y + cos(d);
tmp.mat[2][3] = (1-cos(d))*y*z - x*sin(d);
tmp.mat[3][1] = (1-cos(d))*x*z - y*sin(d);
tmp.mat[3][2] = (1-cos(d))*y*z + x*sin(d);
tmp.mat[3][3] = (1-cos(d))*z*z + cos(d);
tmp.mat[4][4]=1;
p=matMul(tmp,p,4);
return p;
}
Matrix work()
{
Matrix p;
int i;
double x,y,z,d;
memset(p.mat,0,sizeof(p.mat));
for (i=1;i<=4;i++)
p.mat[i][i]=1;
for (i=cnt-3;i>0;i--)
{
if (opstr[i][2]=='T')
{
sscanf(opstr[i],"glTranslatef(%lf,%lf,%lf);",&x,&y,&z);
p=translate(p,x,y,z);
}
else if (opstr[i][2]=='S')
{
sscanf(opstr[i],"glScalef(%lf,%lf,%lf);",&x,&y,&z);
p=scale(p,x,y,z);
}
else if (opstr[i][2]=='R')
{
sscanf(opstr[i],"glRotatef(%lf,%lf,%lf,%lf);",&d,&x,&y,&z);
p=rotate(p,x,y,z,d);
}
}
return p;
}
int main()
{
int t;
double x,y,z;
scanf("%d",&t);
while(t--)
{
cnt=0;
while(1)
{
scanf("%s",opstr[cnt++]);
if (opstr[cnt-1][2]=='E')
break;
}
Matrix ans=work();
sscanf(opstr[cnt-2], "glVertex3f(%lf,%lf,%lf);",&x,&y,&z);
printf("%.1lf %.1lf %.1lf\n",x*ans.mat[1][1]+y*ans.mat[1][2]+z*ans.mat[1][3]+ans.mat[1][4],
x*ans.mat[2][1]+y*ans.mat[2][2]+z*ans.mat[2][3]+ans.mat[2][4],
x*ans.mat[3][1]+y*ans.mat[3][2]+z*ans.mat[3][3]+ans.mat[3][4]);
}
return 0;
}
弄明白了本例，可以順手將HDU 4087 “ALetter to Programmers”做一下，以加深對三維幾何變換的理解。HDU 4087比本例復雜一些，因為存在某個指定的平移、縮放或旋轉操作集重復執行，因此會用到矩陣快速冪運算。
下面給出HDU 4087的源程序供參考。需要說明的是，這個源程序在HDU上用C++提交，一直Wrong Answer。采用G++提交可以Accepted。到現在我也沒弄明白，所以僅供大家參考了。
（3）HDU 4087的參考源程序。
#include <stdio.h>
#include <math.h>
#include <string.h>
#define eps 1e-6
#define pi acos(-1.0)
struct Matrix { double mat[5][5]; // 存儲矩陣中各元素
};
Matrix matMul(Matrix a ,Matrix b,int n) {
Matrix c;
memset(c.mat,0,sizeof(c.mat)); int i,j,k; for (k = 1; k<=n ; k++) for (i=1 ;i<=n ; i++) if (a.mat[i][k]!=0) for (j = 1 ;j<=n ;j++)
c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]); return c; }
Matrix quickMatPow(Matrix a ,int n,int b) // n階矩陣a快速b次冪
{
Matrix c;
memset(c.mat ,0 ,sizeof(c.mat)); int i; for (i = 1 ;i <= n ;i++)
c.mat[i][i] = 1; while (b!=0) { if (b & 1)
c = matMul(c ,a ,n); // c=c*a;
a = matMul(a ,a ,n); // a=a*a
b /= 2; } return c; }
Matrix translate(Matrix p, double x, double y, double z) // 平移
{
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat)); for (int i=1;i<=4;i++)
tmp.mat[i][i]=1;
tmp.mat[1][4] = x;
tmp.mat[2][4] = y;
tmp.mat[3][4] = z;
p=matMul(tmp,p,4); return p; }
Matrix scale(Matrix p, double a, double b, double c) //縮放
{
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat)); for (int i=1;i<=4;i++)
tmp.mat[i][i]=1;
tmp.mat[1][1] = a;
tmp.mat[2][2] = b;
tmp.mat[3][3] = c;
p=matMul(tmp,p,4); return p; }
Matrix rotate(Matrix p, double x, double y, double z, double d) // 旋轉
{
d = d/180*pi; double di = sqrt(x*x+y*y+z*z); //單位化
x/=di, y/=di, z/=di;
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
tmp.mat[1][1] = (1-cos(d))*x*x + cos(d);
tmp.mat[1][2] = (1-cos(d))*x*y - z*sin(d);
tmp.mat[1][3] = (1-cos(d))*x*z + y*sin(d);
tmp.mat[2][1] = (1-cos(d))*x*y + z*sin(d);
tmp.mat[2][2] = (1-cos(d))*y*y + cos(d);
tmp.mat[2][3] = (1-cos(d))*y*z - x*sin(d);
tmp.mat[3][1] = (1-cos(d))*x*z - y*sin(d);
tmp.mat[3][2] = (1-cos(d))*y*z + x*sin(d);
tmp.mat[3][3] = (1-cos(d))*z*z + cos(d);
tmp.mat[4][4]=1;
p=matMul(tmp,p,4); return p; }
Matrix work(int n) {
Matrix p; int i; double x,y,z,d;
memset(p.mat,0,sizeof(p.mat)); for (i=1;i<=4;i++)
p.mat[i][i]=1; char s[100]; while (1) {
scanf("%s", s); if (strcmp(s,"translate")==0) {
scanf("%lf%lf%lf",&x,&y,&z);
p=translate(p,x,y,z); continue; } if (strcmp(s,"scale")==0) {
scanf("%lf%lf%lf",&x,&y,&z);
p=scale(p,x,y,z); continue; } if (strcmp(s,"rotate")==0) {
scanf("%lf%lf%lf%lf",&x,&y,&z,&d);
p=rotate(p,x,y,z,d); continue; } if (strcmp(s,"repeat")==0) { int n1;
scanf("%d", &n1);
p = matMul(work(n1),p,4); } if (strcmp(s,"end")==0) {
p=quickMatPow(p,4,n); return p; } } } int main() { int n,i;
Matrix p; double a[4]; double tx,ty,tz,td; while(scanf("%d", &n) && n!=0) {
p=work(1); while (n--) {
scanf("%lf%lf%lf",&a[1],&a[2],&a[3]);
tx=a[1],ty=a[2],tz=a[3],td=1; for (i=1; i<=4;i++) {
a[i] =p.mat[i][1]*tx +p.mat[i][2]*ty +p.mat[i][3]*tz + p.mat[i][4]*td; if (fabs(a[i])<eps) a[i]=0.0; }
printf("%.2lf %.2lf %.2lf\n",a[1],a[2],a[3]); }
printf("\n"); } return 0; }