矩陣乘法（六）：幾何變換

      在計算機圖形學中，矩陣乘法有着很好的應用。圖形的變換能夠經過構造相應的矩陣進行計算來完成。

      咱們知道，平面上的元素，就是點、線、面，而線就是由一個個點組成的，面是由一條條線組成的，因此歸根結底，平面上全部的圖形都是由點組成的。在座標系中，一個點就是由一對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; }