注意：这是一篇从旧博客恢复的文章。

原地址：http://freemeepo.com/blog/archives/162

---

这是UC_BerkeleyX: CS184.1x 计算机图形学导论的HW1。

因为这两个函数是属于glu而非gl，所以手动实现这两个函数并非没有意义 :-)

考虑到这两个函数关键在于计算出变换矩阵，并将这个矩阵载入到OpenGL中。glLoadMatrixf可以完成载入矩阵的工作。那么矩阵运算我们可以用glm（OpenGL Mathematics）来完成。glm是一个头文件only的库，只需include"glm/glm.hpp"头文件即可。glm可以在这里下载到。需要注意的是，glm里的矩阵是列优先的，即存储方式跟行优先矩阵互为转置。

首先我们写出变换矩阵。

先说gluLookAt的原理：

我们需要从观察者坐标系变换到默认坐标系，需要先平移再旋转。平移矩阵我们很容易写出，只要再左乘一个旋转矩阵即可。由线性代数的知识我们知道，原点不变的坐标系变换矩阵，用该坐标系在原坐标系中的三个正交单位向量组成即可。而且我们知道，用两个向量就可以确定一个坐标系，我们设新坐标系为u, v, w。w就是eye向量，u是up叉乘eye，v是w叉乘u。于是旋转矩阵为：

$$!\begin{bmatrix} x_u & y_u & z_u \\ x_v & y_v & z_v \\ x_w & y_w & z_w \end{bmatrix}$$

但是我们要改成齐次坐标的形式来与平移矩阵相乘，即增加一维，变成：

$$!\begin{bmatrix} x_u & y_u & z_u & 0 \\ x_v & y_v & z_v & 0 \\ x_w & y_w & z_w & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

因此最终结果为：

$$!\begin{bmatrix} x_u & y_u & z_u & 0 \\ x_v & y_v & z_v & 0 \\ x_w & y_w & z_w & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 & -eye_x \\ 0 & 1 & 0 & -eye_y \\ 0 & 0 & 1 & -eye_z \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

不过在我们的这个实现中，center我们定在了原点，并没有考虑center。

gluPerspective的原理：

这个的计算相对复杂一些，中间利用的是相似三角形。具体推导可参考课程。最终可以得到：

$$!\begin{bmatrix} \frac{d}{aspect} & 0 & 0 & 0 \\ 0 & d & 0 & 0 \\ 0 & 0 & A & B \\ 0 & 0 & -1 & 0 \end{bmatrix}$$

其中，$$d=\cot \frac{fovy}{2}$$，$$A=-\frac{zFar+zNear}{zFar-zNear}$$，$$B=-\frac{2\cdot zFar\cdot zNear}{zFar-zNear}$$

代码如下：

#include <windows.h>
#include <GL/glut.h>
#include "glm/glm.hpp"
#include <cstdlib>
#include <cmath>
#include <cstdio>
#define pi (acos(-1.0))
using namespace std;
using namespace glm;

static float c=pi/180.0f;
static int du=90,oldmy=-1,oldmx=-1;
static float r=1.5f,h=0.0f;
mat4 LookAt(const vec3 &eye, const vec3 &up)
{
    mat4 t2(1,0,0,0,0,1,0,0,0,0,1,0,-eye.x,-eye.y,-eye.z,1);
    vec3 w = normalize(eye);
    vec3 u = normalize(cross(up,eye));
    vec3 v = cross(w,u);
    mat4 t1(u.x,v.x,w.x,0,u.y,v.y,w.y,0,u.z,v.z,w.z,0,0,0,0,1);
    return t1*t2;
}
void MygluLookAt(double eyeX,double eyeY,double eyeZ,double centerX,double centerY,double centerZ,double upX,double upY,double upZ)
{
    vec3 eye(eyeX,eyeY,eyeZ);
    vec3 up(upX,upY,upZ);
    mat4 mv=LookAt(eye,up);
    glLoadMatrixf(&mv[0][0]);
}
void display()
{
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    printf("At:%.2f %.2f %.2f\n",r*cos(c*du),h,r*sin(c*du));
    glMatrixMode( GL_MODELVIEW );
//    glLoadIdentity();
//    gluLookAt(r*cos(c*du), h, r*sin(c*du), 0, 0, 0, 0, 1, 0);
    MygluLookAt(r*cos(c*du), h, r*sin(c*du), 0, 0, 0, 0, 1, 0);

    glutWireTeapot(0.5f);

    glutSwapBuffers();
}
void Mouse(int button, int state, int x, int y)
{
    if(state==GLUT_DOWN)
        oldmx=x,oldmy=y;
}
void onMouseMove(int x,int y)
{
    //printf("%d\n",du);
    du+=x-oldmx;
    h +=0.03f*(y-oldmy);
    if(h>1.0f) h=1.0f;
    else if(h<-1.0f) h=-1.0f;
    oldmx=x,oldmy=y;
    glutPostRedisplay();
}
void keyboard(unsigned char key, int x, int y)
{
    if (key==27) exit(0);
}
void MygluPerspective(double fovy, double aspect, double zNear, double zFar)
{
    double d=1/tan(fovy*c/2);
    double A=-(zFar+zNear)/(zFar-zNear);
    double B=-2*zFar*zNear/(zFar-zNear);
    mat4 pers(d/aspect,0,0,0,0,d,0,0,0,0,A,-1,0,0,B,0);
    glLoadMatrixf(&pers[0][0]);
}
void reshape(int w,int h)
{
    glViewport( 0, 0, w, h );
    glMatrixMode( GL_PROJECTION );
//    glLoadIdentity();
//    gluPerspective(75.0f, (float)w/h, 1.0f, 1000.0f);
    MygluPerspective(75.0f, (float)w/h, 1.0f, 1000.0f);
}
int main()
{
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH );
    glutInitWindowSize(800, 600);
    glutCreateWindow("OpenGL");
    glEnable(GL_DEPTH_TEST);
    glutReshapeFunc( reshape );
    glutDisplayFunc(display);
    glutMouseFunc(Mouse);
    glutMotionFunc(onMouseMove);
    glutKeyboardFunc(keyboard);
    glutMainLoop();
    return 0;
}

至于glRotate，实现也不难，只需利用下面的公式：

$$!R(a,\theta)=I_{3\times 3}\cos\theta+aa^T(1-\cos\theta)+A^*\sin\theta$$ $$!R(a,\theta)=\cos\theta\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}+(1-\cos\theta)\begin{bmatrix} x^2 & xy & xz \\ xy & y^2 & yz \\ xz & yz & z^2 \end{bmatrix}+\sin\theta\begin{bmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{bmatrix}$$

需要注意的是，轴a须先单位化。