1. 简介

参考 OpenGL 的逻辑，实现了一个简单的软光栅渲染演示程序，单文件实现，去掉空格和注释 831 行。

模型和纹理无需外部加载，程序启动时自动生成。

模型由正弦函数旋转得到，颜色贴图使用一个改造过的简化版的噪声函数得到，法线贴图和位移贴图从颜色贴图里面计算。

特性：

• 可编程的渲染管线（顶点着色器和片元着色器）
• 环绕相机
• BlinnPhong 着色
• Gouraud 着色
• Flat 着色
• 颜色贴图
• 法线贴图
• 位移贴图
• 线框模式
• 背面剔除 

操作

• 拖动鼠标左键：模型旋转
• 滚轮：视野大小缩放 
• 鼠标左键加滚轮：相机前进或后退（视觉效果和缩放视野差不多）
• 空格键：切换三种光照模型
• 按 C 键：颜色贴图
• 按 N 键：法线贴图
• 按 B 键：位移贴图

2. 效果图

3. 贴图效果

4. 三种光照模型的比对

从左到右依次为 BlinnPhong 着色，Gouraud 着色，Flat 着色。

1.6 万个三角形：BlinnPhong 着色效果最好，但是帧率最低，Flat 着色看起来有点生硬。

9 万个三角形：视觉效果上差别不大，帧率差别不大。

5. 源码实现

/////////////////////////////////////////////////////////////
// Rendering a Vase
// Platform : Visual Studio 2022 (v143), EasyX_20220116
// Author   : 872289455@qq.com
// Date     : 2022-02-21
//
#include <ctime>
#include <cassert>
#include <cfloat>
#include <cmath>
#include <vector>
#include <random>
#include <functional>
#include <algorithm>
#include <easyx.h>
using std::vector;

#define WIDTH                800
#define HEIGHT               600
#define Z_NEAR               (1.f)
#define Z_FAR                (50.f)
#define WORLD_UP             vec3(0.f, 1.f, 0.f)
#define MODE_FILL            0
#define MODE_LINE            1
#define COLOR_STRENGTH       0.85f
#define NORMAL_STRENGTH      3.f
#define DISPACE_STRENGTH     0.23f
#define RADIANS(degrees)     ((degrees) * 0.0174532925f)
#define DEGREES(radians)     ((radians) * 57.295779513f)
#define M_PI                 3.14159265358979323846f

///////////////////////////////////////////////////////////////////////////////////
// Math
///////////////////////////////////////////////////////////////////////////////////
struct vec2 {
    float x, y;
    vec2(float a = 0, float b =0) : x(a), y(b) {}
    vec2 operator*(float f) {
        return vec2(f * x, f * y);
    }
    vec2 operator-(vec2 b) {
        return vec2(x-b.x, y-b.y);
    }
    vec2 operator+(vec2 b) {
        return vec2(x + b.x, y + b.y);
    }
};

struct vec3 {
    float x, y, z;
    vec3(float a = 0, float b = 0, float c = 0) : x(a), y(b), z(c){}
    vec3 operator*(float f) {
        return vec3(f * x, f * y, f * z);
    }
    vec3 operator-() const {
        return vec3(-x, -y, -z);
    }
    vec3 operator-(float f) const {
        return vec3(x - f, y - f, z - f);
    }
    vec3 operator-(vec3 b) const {
        return vec3(x - b.x, y - b.y, z - b.z);
    }
    vec3 operator+(vec3 b) const {
        return vec3(x + b.x, y + b.y, z + b.z);
    }
    float length() {return sqrtf(x*x+y*y+z*z);}
};

struct vec4 {
    float x, y, z, w;
    vec4(float a = 0, float b = 0, float c = 0, float d = 0) : x(a), y(b), z(c), w(d) {}
    vec4(vec3 v, float f = 1.f) : x(v.x), y(v.y), z(v.z), w(f) {}
    vec4 operator*(float f) {
        return vec4(f * x, f * y, f * z, f * w);
    }
    vec3 to_vec3() {
        return vec3(x, y, z);
    }
};

vec2 operator*(float f, vec2 vec) {
    return vec * f;;
}

vec3 operator*(float f, vec3 vec) {
    return vec * f;;
}

vec4 operator*(float f, vec4& vec) {
    return vec * f;;
}

float dot(vec3 a, vec3 b) {
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

vec3 cross(vec3 a, vec3 b) {
    return vec3(a.y * b.z - b.y * a.z,
                a.z * b.x - b.z * a.x,
                a.x * b.y - b.x * a.y
                );
}

vec2 normalize(vec2 v) {
    float len_sq = v.x * v.x + v.y * v.y;
    if (len_sq < FLT_EPSILON)
        return v;
    return (1.f / sqrt(len_sq)) * v;
}

vec3 normalize(vec3 v) {
    float len_sq = v.x * v.x + v.y * v.y + v.z * v.z;
    if (len_sq < FLT_EPSILON)
        return v;
    return (1.f / sqrt(len_sq)) * v;
}

struct quat {
    union {
        struct { float x, y, z, w; };
        struct { vec3 vec; float scalar; };
    };
    quat(float angle, vec3 &axis) {
        float s = 0.f;
        float c = 0.f;
        float rads = RADIANS(angle);
        axis = normalize(axis);
        s = sinf(rads * 0.5f);
        c = cosf(rads * 0.5f);
        x = s * axis.x;
        y = s * axis.y;
        z = s * axis.z;
        w = c;
    }
    vec3 rotate(vec3 from) {
        vec3 a = 2.f * dot(vec, from) * vec;
        vec3 b = (w * w - dot(vec, vec)) * from;
        vec3 c = 2.f * w * cross(vec, from);
        return a + b + c;
    }
};

struct mat3 {
    vec3 col[3];
    mat3(vec3 c0, vec3 c1, vec3 c2) {
        col[0] = c0;
        col[1] = c1;
        col[2] = c2;
    }
    vec3 operator*(vec3 v) {
        return vec3(
            col[0].x * v.x + col[1].x * v.y + col[2].x * v.z,
            col[0].y * v.x + col[1].y * v.y + col[2].y * v.z,
            col[0].z * v.x + col[1].z * v.y + col[2].z * v.z
        );
    }
};

struct mat4 {
    union {
        float v[16];
        struct { vec4 col[4]; };
    };
    mat4(vec4 c0 = vec4(1.f, 0.f, 0.f, 0.f),
         vec4 c1 = vec4(0.f, 1.f, 0.f, 0.f),
         vec4 c2 = vec4(0.f, 0.f, 1.f, 0.f),
         vec4 c3 = vec4(0.f, 0.f, 0.f, 1.f)
        ){
        col[0] = c0;
        col[1] = c1;
        col[2] = c2;
        col[3] = c3;
    }
    vec4& operator[](int i) { return col[i]; }
    mat4  operator*(float f) {
        return mat4(f * col[0], f * col[1], f * col[2], f * col[3]);
    }
    vec4 operator*(vec4 v) {
        return vec4(
            col[0].x * v.x + col[1].x * v.y + col[2].x * v.z + col[3].x * v.w,
            col[0].y * v.x + col[1].y * v.y + col[2].y * v.z + col[3].y * v.w,
            col[0].z * v.x + col[1].z * v.y + col[2].z * v.z + col[3].z * v.w,
            col[0].w * v.x + col[1].w * v.y + col[2].w * v.z + col[3].w * v.w
        );
    }
    mat4 operator*(mat4 &m) {
        return mat4(
            *this * m.col[0],
            *this * m.col[1],
            *this * m.col[2],
            *this * m.col[3]
        );
    }
};

mat4 transposed(mat4 m) {
    return mat4(vec4(m[0].x, m[1].x, m[2].x, m[3].x),
                vec4(m[0].y, m[1].y, m[2].y, m[3].y),
                vec4(m[0].z, m[1].z, m[2].z, m[3].z),
                vec4(m[0].w, m[1].w, m[2].w, m[3].w)
                );
}

#define M4_3X3MINOR(c0, c1, c2, r0, r1, r2) \
    (m.v[c0 * 4 + r0] * (m.v[c1 * 4 + r1] * m.v[c2 * 4 + r2] - m.v[c1 * 4 + r2] * m.v[c2 * 4 + r1]) - \
     m.v[c1 * 4 + r0] * (m.v[c0 * 4 + r1] * m.v[c2 * 4 + r2] - m.v[c0 * 4 + r2] * m.v[c2 * 4 + r1]) + \
     m.v[c2 * 4 + r0] * (m.v[c0 * 4 + r1] * m.v[c1 * 4 + r2] - m.v[c0 * 4 + r2] * m.v[c1 * 4 + r1]))

float determinant(mat4 &m) {
    return  m.v[0] * M4_3X3MINOR(1, 2, 3, 1, 2, 3)
          - m.v[4] * M4_3X3MINOR(0, 2, 3, 1, 2, 3)
          + m.v[8] * M4_3X3MINOR(0, 1, 3, 1, 2, 3)
          - m.v[12]* M4_3X3MINOR(0, 1, 2, 1, 2, 3);
}

mat4 adjugate(mat4 &m) {
    mat4 ret;
    float* cv = &ret.v[0];
    cv[0] = M4_3X3MINOR(1, 2, 3, 1, 2, 3);
    cv[1] =-M4_3X3MINOR(1, 2, 3, 0, 2, 3);
    cv[2] = M4_3X3MINOR(1, 2, 3, 0, 1, 3);
    cv[3] =-M4_3X3MINOR(1, 2, 3, 0, 1, 2);
    cv[4] =-M4_3X3MINOR(0, 2, 3, 1, 2, 3);
    cv[5] = M4_3X3MINOR(0, 2, 3, 0, 2, 3);
    cv[6] =-M4_3X3MINOR(0, 2, 3, 0, 1, 3);
    cv[7] = M4_3X3MINOR(0, 2, 3, 0, 1, 2);
    cv[8] = M4_3X3MINOR(0, 1, 3, 1, 2, 3);
    cv[9] =-M4_3X3MINOR(0, 1, 3, 0, 2, 3);
    cv[10]= M4_3X3MINOR(0, 1, 3, 0, 1, 3);
    cv[11]=-M4_3X3MINOR(0, 1, 3, 0, 1, 2);
    cv[12]=-M4_3X3MINOR(0, 1, 2, 1, 2, 3);
    cv[13]= M4_3X3MINOR(0, 1, 2, 0, 2, 3);
    cv[14]=-M4_3X3MINOR(0, 1, 2, 0, 1, 3);
    cv[15]= M4_3X3MINOR(0, 1, 2, 0, 1, 2);
    return transposed(ret);
}

mat4 inverse(mat4 m) {
    float det = determinant(m);
    if (fabs(det) < FLT_EPSILON)
        return mat4();
    return adjugate(m) * (1.f / det);
}

///////////////////////////////////////////////////////////////////////////////////
// Modeling a vase : x = 2 + sin(y)
///////////////////////////////////////////////////////////////////////////////////
struct Vase {
    vector<vec3> pos;
    vector<vec3> normal;
    vector<vec3> tangent;
    vector<vec2> uv;
    vector<int> index;

    Vase(int y_edges) {
        int x_edges = 2 * y_edges;
        float dx = 2.f * M_PI / x_edges;
        float dy = 2.f * M_PI / y_edges;
        for (int y = 0; y <= y_edges; y++) {
            float v = y * dy;
            float cv = cos(v);
            float sv = sin(v);
            float r = 2 + sv;
            for (int x = 0; x <= x_edges; x++) {
                float u = x * dx;
                float cu = cos(u);
                float su = sin(u);
                vec3 tang(cu * cv, 1, su * cv);
                vec3 biTangent(-su, 0, cu);
                pos.push_back(vec3(r * cu, v, r * su));
                tangent.push_back(tang);
                normal.push_back(normalize(cross(tang, biTangent)));
                uv.push_back(vec2(u / (2 * M_PI), v / (2 * M_PI)));
            }
        }
        for (int y = 0; y < y_edges; y++) {
            int start = y * (x_edges + 1);
            int next_start = (y + 1) * (x_edges + 1);
            for (int x = 0; x < x_edges; x++) {
                int next_x = x + 1;
                index.push_back(start + x);
                index.push_back(next_start + x);
                index.push_back(start + next_x);
                index.push_back(start + next_x);
                index.push_back(next_start + x);
                index.push_back(next_start + next_x);
            }
        }
        pos.push_back(vec3());
        tangent.push_back(vec3(1.f, 0.f, 0.f));
        normal.push_back(vec3(0.f, -1.f, 0.f));
        uv.push_back(vec2(0.5, 0.5));
        for (int x = 0; x < x_edges; x++) {
            index.push_back((int)pos.size() - 1);
            index.push_back(x);
            index.push_back(x + 1);
            normal[x] = vec3(0.f, -1.f, 0.f);
        }
        normal[x_edges] = vec3(0.f, -1.f, 0.f);
    }
};

///////////////////////////////////////////////////////////////////////////////////
// Procedural Texturing
///////////////////////////////////////////////////////////////////////////////////
class Noise {
    const uint32_t TABLE_SIZE = 0X100U;
    const uint32_t TABLE_MASK = 0XFFU;
    vector<vector<vec2>> grid;
public:
    Noise() {
        std::mt19937 gen(time(nullptr));
        std::uniform_real_distribution<float> distrFloat(-1.f, 1.f);
        auto randFloat = std::bind(distrFloat, gen);
        grid = vector<vector<vec2>> (TABLE_SIZE, vector<vec2>(TABLE_SIZE));
        for (uint32_t row = 0; row < TABLE_SIZE; row++)
            for (uint32_t col = 0; col < TABLE_SIZE; col++) 
                grid[row][col] = normalize(vec2(randFloat(), randFloat()));
    }

    vec3 eval(vec2 p) {
        uint32_t x = (uint32_t)std::floor(p.x);
        uint32_t y = (uint32_t)std::floor(p.y);
        float u = fade(p.x - x);
        float v = fade(p.y - y);
        uint32_t x0 = x & TABLE_MASK;
        uint32_t x1 = (x0 + 1) & TABLE_MASK;
        uint32_t y0 = y & TABLE_MASK;
        uint32_t y1 = (y0 + 1) & TABLE_MASK;
        vec2 vec = lerp(lerp(grid[x0][y0], grid[x1][y0], u), lerp(grid[x0][y1], grid[x1][y1], u), v);
        return vec3(vec.x, vec.y, COLOR_STRENGTH);
    }
private:
    vec2 lerp(vec2 lo, vec2 hi, float t) {
        return lo * (1 - t) + hi * t;
    }
    float fade(float t) {
        return t * t * t * (t * (t * 6 - 15) + 10);
    }
};

struct Image {
    int height;
    int width;
    int channels;
    vector<uint8_t> data;
    Image(int h = 0, int w = 0, int c = 0, vector<uint8_t> raw = vector<uint8_t>()) : height(h), width(w), channels(c), data(raw) {}
};

vec3 getPixelAt(const Image& map, int x, int y) {
    if (x < 0) x = map.width - 1;
    if (x >= map.width) x = 0;
    if (y < 0) y = map.height - 1;;
    if (y >= map.height) y = 0;

    float r, g, b;
    int ofs = map.channels * (y * map.width + x);
    r = g = b = map.channels > 0 ? map.data[ofs + 0] : 0.f;
    g = b = map.channels > 1 ? map.data[ofs + 1] : r;
    b = map.channels > 2 ? map.data[ofs + 2] : g;
    return (1.f / 255.f) * vec3(r, g, b);
}

// vec3 returned in range 0.0 - 1.0
vec3 texture(const Image& map, const vec2& texCoord) {
    vec2 uv = texCoord;
    int x = int(map.width * uv.x);
    int y = int(map.height * uv.y);
    float u = map.width * uv.x - x;
    float v = map.height * uv.y - y;

    vec3 c00 = getPixelAt(map, x, y);
    vec3 c10 = getPixelAt(map, x + 1, y);
    vec3 c01 = getPixelAt(map, x, y + 1);
    vec3 c11 = getPixelAt(map, x + 1, y + 1);
    return (1 - u) * (1 - v) * c00 + u * (1 - v) * c10 + (1 - u) * v * c01 + u * v * c11;
}

Image buildColorTexture(void) {
    Noise noise;
    float baseFreq = 12;
    int level = 6;
    float amplitudeMult = 0.5f;
    float frequencyMult = 1.9f;
    int width = 1024;
    int height = 1024;

    vector<uint8_t> pixels;
    pixels.resize(width * height * 3);

    float dx = 1.0f / (width - 1);
    float dy = 1.0f / (height - 1);
    for (int row = 0; row < height; row++) {
        for (int col = 0; col < width / 2; col++) {
            float x = dx * col;
            float y = dy * row;
            float freq = baseFreq;
            float amplitude = 1.f;
            vec3 sumVec;
            for (int oct = 0; oct < level; oct++) {
                sumVec = sumVec + amplitude * noise.eval(freq * vec2(x, y));
                freq *= frequencyMult;
                amplitude *= amplitudeMult;
            }
            sumVec = normalize(sumVec);
            float turb = sin(sumVec.y);
            int lofs = 3 * (row * width + col);
            int rofs = 3 * (row * width + width - 1 - col);
            pixels[lofs + 0] = pixels[rofs + 0] = static_cast<uint8_t>((turb * 1.0 + 1.0) * 255.f * 0.5);
            pixels[lofs + 1] = pixels[rofs + 1] = static_cast<uint8_t>((turb * 0.7 + 1.0) * 255.f * 0.5);
            pixels[lofs + 2] = pixels[rofs + 2] = static_cast<uint8_t>((turb * 0.8 + 1.0) * 255.f * 0.5);
        }
    }
    return Image(width, height, 3, pixels);
}

Image buildNormalTextureFrom(const Image& src) {
    vector<uint8_t> pixels;
    pixels.reserve(3 * src.width * src.height);
    for (int y = 0; y < src.height; y++) {
        for (int x = 0; x < src.width; x++) {
            float tl = getPixelAt(src, x - 1, y - 1).length();
            float t = getPixelAt(src, x - 1, y).length();
            float tr = getPixelAt(src, x - 1, y + 1).length();
            float r = getPixelAt(src, x, y + 1).length();
            float br = getPixelAt(src, x + 1, y + 1).length();
            float b = getPixelAt(src, x + 1, y).length();
            float bl = getPixelAt(src, x + 1, y - 1).length();
            float l = getPixelAt(src, x, y - 1).length();
            float dx = float(tr + 2.f * r + br) - float(tl + 2.f * l + bl);
            float dy = float(bl + 2.f * b + br) - float(tl + 2.f * t + tr);
            float dz = 1.f / NORMAL_STRENGTH;
            vec3 n = normalize(vec3(dx, dy, dz));
            pixels.push_back(static_cast<uint8_t>((n.x + 1.0) * 255.f * 0.5));
            pixels.push_back(static_cast<uint8_t>((n.y + 1.0) * 255.f * 0.5));
            pixels.push_back(static_cast<uint8_t>((n.z + 1.0) * 255.f * 0.5));
        }
    }
    return Image(src.width, src.height, 3, pixels);
}


///////////////////////////////////////////////////////////////////////////////////
// Orbit Camera
///////////////////////////////////////////////////////////////////////////////////
uint8_t activeColorTexture = 1u;
uint8_t activeNormalTexture = 0u;
uint8_t activeDispTexture = 0u;
uint8_t activeLineMode = 0u;
uint8_t cullFace = 0u;
uint8_t selectShader = 0u;

class Camera {
    vec3  pos, target;
    float zoom, aspect;
    mat4  projection, view;
public:
    Camera(vec3 position, vec3 center, float asp) : pos(position), target(center), zoom(25.f), aspect(asp) {
        updateProjection();
        updateView();
    }
    void processInput() {
        ExMessage msg;
        static bool firstClick = true;
        static int lastx, lasty;
        while (peekmessage(&msg, EM_MOUSE | EM_KEY)) {
            // Drag Left Button: Rotate
            if (WM_MOUSEMOVE == msg.message) {
                if (false == msg.lbutton) {
                    firstClick = true;
                }
                else if (firstClick) {
                    firstClick = false;
                    lastx = msg.x;
                    lasty = msg.y;
                }
                else {
                    float xoff = 0.5f * (msg.x - lastx);
                    float yoff = 0.5f * (msg.y - lasty);
                    lastx = msg.x;
                    lasty = msg.y;
                    mouseMove(xoff, yoff);
                }
            }
            else if (WM_MOUSEWHEEL == msg.message) {
                if (msg.lbutton)    // Left Button + Scroll: Camera Forward or Backward, have a bug
                    srcollPosition(-0.01f * msg.wheel);
                else                // Scroll: field of view
                    scrollZoom(0.01f * msg.wheel);
            }
            else if (WM_KEYDOWN == msg.message) {
                switch (msg.vkcode) {
                case VK_SPACE: selectShader = (selectShader + 1) % 3;    break;
                case 'C':      activeColorTexture ^= 1u;                 break;
                case 'N':      activeNormalTexture ^= 1u;                break;
                case 'B':      activeDispTexture ^= 1u;                  break;
                case 'L':      activeLineMode ^= 1u;                     break;
                case 'K':      cullFace ^= 1u;                           break;
                default:                                                 break;
                }
            }
        }
    }
    mat4& viewMatrix()        { return view; }
    mat4& projectionMatrix()  { return projection;}
    vec3& position()          { return pos; }
    float fov()               { return zoom; }
private:
    void updateProjection() {
        float t = Z_NEAR * tan(RADIANS(0.5f * zoom));
        float b = -t;
        float r = t * aspect;
        float l = -r;
        projection[0] = vec4(2 * Z_NEAR / (r - l), 0.f, 0.f, 0.f);
        projection[1] = vec4(0.f, 2 * Z_NEAR / (t - b), 0.f, 0.f);
        projection[2] = vec4((r + l) / (r - l), (t + b) / (t - b), -(Z_FAR - Z_NEAR) / (Z_FAR - Z_NEAR), -1.f);
        projection[3] = vec4(0.f, 0.f, -2 * Z_FAR * Z_NEAR / (Z_FAR - Z_NEAR), 0.f);
    }
    void updateView() {
        vec3 z_axis = normalize(pos - target);
        vec3 x_axis = normalize(cross(WORLD_UP, z_axis));
        vec3 y_axis = normalize(cross(z_axis, x_axis));

        view[0] = vec4(x_axis.x, y_axis.x, z_axis.x, 0.f);
        view[1] = vec4(x_axis.y, y_axis.y, z_axis.y, 0.f);
        view[2] = vec4(x_axis.z, y_axis.z, z_axis.z, 0.f);
        view[3] = vec4(-dot(x_axis, pos), -dot(y_axis, pos), -dot(z_axis, pos), 1.f);
    }

    void mouseMove(float xoff, float yoff) {
        vec3 front = normalize(target - pos);
        vec3 right = normalize(cross(front, WORLD_UP));
        vec3 up = normalize(cross(right, front));
        float angle = DEGREES(acos(dot(front, WORLD_UP)));
        if ((yoff < 0 && (angle + yoff < 20)) || (yoff > 0 && (angle + yoff > 160)))
            return;

        quat qx(-xoff, up);
        pos = qx.rotate(pos - target) + target;
        front = normalize(target - pos);
        right = normalize(cross(front, up));

        quat qy(-yoff, right);
        pos = qy.rotate(pos - target) + target;
        updateView();
    }

    void srcollPosition(float yoff) {
        vec3 front = normalize(target - pos);
        pos = pos + yoff * front;
        updateView();
    }
    
    void scrollZoom(float yoff) {
        if (zoom + yoff > 120 || zoom + yoff < 10)
            return;
        zoom += yoff;
        updateProjection();
    }
};

///////////////////////////////////////////////////////////////////////////////////
// Rendering
///////////////////////////////////////////////////////////////////////////////////
struct Vertex {
    vec3 pos, normal, tangent, color;
    vec2 uv;
};

typedef void (*VERTEX_SHADER)(const Vertex& in, Vertex& out, vec4& mvpPos);
typedef vec3(*FRAGMENT_SHADER)(const Vertex& in);

struct Program {
    bool              flatHint;
    VERTEX_SHADER     vertexShader;
    FRAGMENT_SHADER   fragmentShader;
    Program(VERTEX_SHADER vert = nullptr, FRAGMENT_SHADER frag = nullptr, bool isFlat = false) : vertexShader(vert), fragmentShader(frag), flatHint(isFlat) {}
};

struct Light {
    vec3  dir;      // Light direction in eye coords.
    float La;       // Ambient light intensity
    float Ld;       // Diffuse light intensity
    float Ls;       // Specular light intensity
    Light(vec3 d = vec3(0.f, 1.f, 1.f), float la = 0.03f, float ld = 0.2f, float ls = 0.8f) : dir(d), La(la), Ld(ld), Ls(ls) {}
};

// Sharing of uniform data between shader programs.
struct perFrameData {
    mat4 modelViewMatrix;
    mat4 projectionMatrix;
    mat4 normalMatrix;
    mat4 viewportMatrix;
    mat4 MVP;

    Image colorMap;
    Image normalMap;
    Image dispMap;

    Light light[4] = {
        Light(vec3( 1.0f,  1.f, 1.f)), 
        Light(vec3(-2.0f, -2.f, 1.f)), 
        Light(vec3(-0.5f,  0.f, 1.f)), 
        Light(vec3(-1.0f,  0.f, 0.f)) 
    };
} uniform;

void setMatrixUniform(mat4 &v, mat4 &proj) {
    float w2 = WIDTH / 2.f;
    float h2 = HEIGHT / 2.f;
    uniform.modelViewMatrix = v;
    uniform.projectionMatrix = proj;
    uniform.normalMatrix = transposed(inverse(v));
    uniform.viewportMatrix = mat4(vec4(w2, 0.f, 0.f, 0.f), vec4(0.f, h2, 0.f, 0.f), vec4(0.f, 0.f, 1.f, 0.f), vec4(w2 - 0.5f, h2 - 0.5f, 0.f, 1.f));
    uniform.MVP = proj * v;
}

class Render {
    int             vertices;
    DWORD*          framebuf;
    VERTEX_SHADER   vertexShader;
    FRAGMENT_SHADER fragmentShader;
    vector<Vertex>  VertexIn;
    vector<Vertex>  VertexOut;
    vector<vec4>    mvpPos;
    vector<int>     index;
    vector<float>   zbuf;
    int             drawMode;
    bool            flatHint;

public:
    Render() : vertices(0), vertexShader(nullptr), fragmentShader(nullptr), framebuf(GetImageBuffer()), drawMode(MODE_FILL), flatHint(false) { assert(framebuf); }

    void createBuffersWith(Vase *vase) {
        assert(vase);
        index = vase->index;
        vertices = (int)vase->pos.size();
        mvpPos.resize(vertices);
        VertexIn.resize(vertices);
        VertexOut.resize(vertices);
        zbuf.resize(WIDTH * HEIGHT);
        for (int i = 0; i < vertices; i++) {
            VertexIn[i].pos = vase->pos[i];
            VertexIn[i].normal = vase->normal[i];
            VertexIn[i].tangent = vase->tangent[i];
            VertexIn[i].uv = vase->uv[i];
        }
    }

    void useProgram(Program& shader) {
        flatHint = shader.flatHint;
        vertexShader = shader.vertexShader;
        fragmentShader = shader.fragmentShader;
    }

    void drawElements(void) {
        std::fill_n(zbuf.begin(), zbuf.size(), FLT_MAX);
        #pragma omp parallel for
        // 1. Vertex Shader
        for (int i = 0; i < vertices; i++)
            vertexShader(VertexIn[i], VertexOut[i], mvpPos[i]);
        
        #pragma omp parallel for
        for (int i = 0; i < (int)index.size(); i += 3) {
            // 2. Shape Assembly
            Vertex vert[3] = { VertexOut[index[i]], VertexOut[index[i + 1]], VertexOut[index[i + 2]] };
            vec4 v[3] = {mvpPos[index[i]], mvpPos[index[i + 1]], mvpPos[index[i + 2]] };
            // 3. Rasterization
            for (auto& p : v) {    
                p = uniform.viewportMatrix * p;
                p = p * (1.f / p.w);
            }
            
            if (cullFace && cross(v[1].to_vec3() - v[0].to_vec3(), v[2].to_vec3() - v[0].to_vec3()).z <= 0) 
                continue;
            
            if (activeLineMode) {
                drawLine(v[0], v[1]);
                drawLine(v[0], v[2]);
                drawLine(v[1], v[2]);
                continue;
            }
            
            int left = (int)min(v[0].x, min(v[1].x, v[2].x)) - 1;
            int right = (int)max(v[0].x, max(v[1].x, v[2].x)) + 1;
            int bottom = (int)min(v[0].y, min(v[1].y, v[2].y)) - 1;
            int top = (int)max(v[0].y, max(v[1].y, v[2].y)) + 1;

            for (int x = left; x <= right; x++) {
                if (x >= WIDTH ||  x <= 0) 
                    continue;
                for (int y = bottom; y <= top; y++) {
                    if (y >= HEIGHT || y <= 0)
                        continue;
                    vec3 centric = barycentric(x, y, v);
                    if (insideTriangle(centric)) {
                        if (flatHint) {
                            float z = v[1].z;
                            if (z < zbuf[indexOf(x, y)]) {
                                zbuf[indexOf(x, y)] = z;
                                // 4. Fragment Shader, flat shading use the second point of triangle
                                setpixel(x, y, fragmentShader(vert[1]));
                            }
                        }
                        else {
                            float a = centric.x, b = centric.y, c = centric.z;
                            float z = a * v[0].z + b * v[1].z + c * v[2].z;
                            if (z < zbuf[indexOf(x, y)]) {
                                Vertex fragIn;
                                zbuf[indexOf(x, y)] = z;
                                fragIn.pos = a * vert[0].pos + b * vert[1].pos + c * vert[2].pos;
                                fragIn.normal = a * vert[0].normal + b * vert[1].normal + c * vert[2].normal;
                                fragIn.tangent = a * vert[0].tangent + b * vert[1].tangent + c * vert[2].tangent;
                                fragIn.color = a * vert[0].color + b * vert[1].color + c * vert[2].color;
                                fragIn.uv = a * vert[0].uv + b * vert[1].uv + c * vert[2].uv;
                                // 4. Fragment Shader
                                setpixel(x, y, fragmentShader(fragIn));
                            }
                        }

                    }
                }
            }

        }
    }
    void bindTexture(const Image &colorMap, const Image &normalMap, const Image & dispMap) {
        uniform.colorMap = colorMap;
        uniform.normalMap= normalMap;
        uniform.dispMap  = dispMap;
    }
    void polygonMode(int drawhint) {
        drawMode = drawhint;
    }

private:
    int indexOf(int x, int y) {
        return (HEIGHT - y) * WIDTH + x;
    }

    void setpixel(int x, int y, vec3 fragColor) {
        uint8_t r = (uint8_t)max(0.f, min(255.f, 255.f * fragColor.x));
        uint8_t g = (uint8_t)max(0.f, min(255.f, 255.f * fragColor.y));
        uint8_t b = (uint8_t)max(0.f, min(255.f, 255.f * fragColor.z));
        framebuf[indexOf(x, y)] = 0Xffu << 24 | r << 16 | g << 8 | b; // BGRA
    }

    bool insideTriangle(vec3& centric) {
        return  centric.x >= 0.f && centric.x <= 1.f && \
                centric.y >= 0.f && centric.y <= 1.f && \
                centric.z >= 0.f && centric.z <= 1.f;
    }
     
    vec3 barycentric(int x, int y, const vec4* _v) {
        vec2 p(static_cast<float>(x), static_cast<float>(y));
        vec2 a(_v[0].x, _v[0].y), b(_v[1].x, _v[1].y), c(_v[2].x, _v[2].y);
        vec2 v0 = b - a, v1 = c - a, v2 = p - a;
        float div = 1.f / (v0.x * v1.y - v1.x * v0.y);
        float v = div * (v2.x * v1.y - v1.x * v2.y);
        float w = div * (v0.x * v2.y - v2.x * v0.y);
        float u = 1.f - v - w;
        return vec3(u,v,w);
    }

    void drawLine(vec4 begin, vec4 end) {
        vec3 lineColor(0.45F, 0.73F, 0.19F);
        int x0 = int(begin.x);
        int y0 = int(begin.y);
        int x1 = int(end.x);
        int y1 = int(end.y);
        if (x0 < 0 || x0 >= WIDTH || y0 < 0 || y0 >= HEIGHT || x1 < 0 || x1 >= WIDTH || y1 < 0 || y1 >= HEIGHT)
            return;
        int dx = abs(x1 - x0), sx = x0 < x1 ? 1 : -1;
        int dy = abs(y1 - y0), sy = y0 < y1 ? 1 : -1;
        int err = (dx > dy ? dx : -dy) / 2;
        while (setpixel(x0, y0, lineColor), x0 != x1 || y0 != y1) {
            int e2 = err;
            if (e2 > -dx) { err -= dy; x0 += sx; }
            if (e2 < dy) { err += dx; y0 += sy; }
        }
    }
};


///////////////////////////////////////////////////////////////////////////////////
// Shaders
///////////////////////////////////////////////////////////////////////////////////
// Common functions for all shaders
vec3 updatePosition(const vec3& old, const vec3& normal, const vec2& uv) {
    vec3 out = old;
    vec3 h = texture(uniform.dispMap, uv);
    float hs = (h.x + h.y + h.z) / 3.f;
    float t = 2 * hs - 1;
    return old - DISPACE_STRENGTH * t * normal;
}

vec3 updateNormal(vec3 Normal, vec3 Tangent, vec2 uv) {
    vec3 N = normalize(Normal);
    vec3 T = normalize(Tangent - N * dot(Tangent, N));
    vec3 B = normalize(cross(N, T)); 
    mat3 TBN(T, B, N);
    vec3 normal = 2.0 * texture(uniform.normalMap, uv) - 1.f;
    return normalize(TBN * normal);
}

vec3 blinnPhongColor(const vec3 &pos, const vec3 &n, const vec2 &uv) {
    vec3 color(0.3f, 0.3f, 0.3f);
    if (activeColorTexture)
        color = texture(uniform.colorMap, uv);
    float la = 0.f;
    vec3 ambient, diffuse, specular;
    for (auto &l : uniform.light) {
        ambient = ambient + l.La * color;
        vec3 s = l.dir;
        float sDotN = max(dot(s, n), 0.f);
        diffuse = diffuse + l.Ld * sDotN * color;
        if (sDotN > 0.0) {
            vec3 v = normalize(-pos);
            vec3 h = normalize(v + s);
            specular = specular + l.Ls * powf(max(dot(h, n), 0), 128) * vec3(0.8f, 0.8f, 0.8f);
        }
    }
    return ambient + diffuse + specular;
}

///////////////////////////////////////////////////////////////////////////////////
// 1. Blinn-Phong Shading: shade each pixel
///////////////////////////////////////////////////////////////////////////////////
void phongVertShader(const Vertex& in, Vertex& out, vec4& mvpPos) {
    vec3 pos = in.pos;
    if (activeDispTexture)
        pos = updatePosition(in.pos, in.normal, in.uv);
    out.pos = (uniform.modelViewMatrix * vec4(pos)).to_vec3();
    out.normal = normalize((uniform.normalMatrix * vec4(in.normal, 0.f)).to_vec3());
    out.tangent = normalize((uniform.normalMatrix * vec4(in.tangent, 0.f)).to_vec3());
    out.uv = in.uv;
    mvpPos = uniform.MVP * vec4(pos);
}

vec3 phongFragShader(const Vertex& in) {
    vec3 N = normalize(in.normal);
    if (activeNormalTexture)
        N = updateNormal(N, in.tangent, in.uv);
    // two-sided shading
    vec3 view = normalize(-in.pos);
    float vDotN = dot(view, N);
    if (vDotN <= 0)
       N = -N;
    return blinnPhongColor(in.pos, N, in.uv);
}


///////////////////////////////////////////////////////////////////////////////////
// 2. Gouraud Shading: shade each vertext
///////////////////////////////////////////////////////////////////////////////////
void gouraudVertShader(const Vertex& in, Vertex& out, vec4& mvpPos) {
    vec3 pos = in.pos;
    if (activeDispTexture)
        pos = updatePosition(in.pos, in.normal, in.uv);
    pos = (uniform.modelViewMatrix * vec4(pos)).to_vec3();
    vec3 N = normalize((uniform.normalMatrix * vec4(in.normal, 0.f)).to_vec3());
    vec3 T = normalize((uniform.normalMatrix * vec4(in.tangent, 0.f)).to_vec3());
    if (activeNormalTexture)
        N = updateNormal(N, T, in.uv);
    vec3 view = normalize(-pos);
    float vDotN = dot(view, N);
    if (vDotN < 0)
       N = -N;
    out.color = blinnPhongColor(pos, N, in.uv);
    mvpPos = uniform.projectionMatrix * vec4(pos);
}

vec3 gouraudFragShader(const Vertex& in) {
    return in.color;
}

///////////////////////////////////////////////////////////////////////////////////
// 3. Flat Shading: shade each triangle
///////////////////////////////////////////////////////////////////////////////////
void flatVertShader(const Vertex& in, Vertex& out, vec4& mvpPos) {
    return gouraudVertShader(in, out, mvpPos);
}
vec3 flatFragShader(const Vertex& in) {
    return gouraudFragShader(in);
}

int main() {
    TCHAR  buf[128];
    size_t vertices, triangles;
    double deltaTime = 0.f;
    double lastTime = 0.f;
    double currentTime = 0.f;
    uint64_t start, now, frequency;

    initgraph(WIDTH, HEIGHT);

    Render render;
    Camera camera(vec3(0.f, 9.f, 19.f), vec3(0.f, M_PI, 0.f), (float)WIDTH / (float)HEIGHT);

    Vase *vase = new Vase(64);
    render.createBuffersWith(vase);
    vertices = vase->pos.size();
    triangles = vase->index.size() / 3;
    delete vase;

    Program blinnPhong(phongVertShader, phongFragShader);
    Program gouraud(gouraudVertShader, gouraudFragShader);
    Program flat(flatVertShader, flatFragShader, true);
    Program shaders[3] = {blinnPhong, gouraud, flat};
    const TCHAR  *shaderName[3] = {_T("blinnPhong"), _T("gouraud"), _T("flat")};

    Image colorMap = buildColorTexture();
    Image normalMap = buildNormalTextureFrom(colorMap);
    render.bindTexture(colorMap, normalMap, colorMap);

    QueryPerformanceCounter((LARGE_INTEGER*)&start);
    QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
    BeginBatchDraw();
    while (1) {
        QueryPerformanceCounter((LARGE_INTEGER*)&now);
        currentTime = (double)(now - start) / frequency;
        deltaTime = currentTime - lastTime;
        lastTime = currentTime;

        camera.processInput();
        setMatrixUniform(camera.viewMatrix(), camera.projectionMatrix());

        if (activeLineMode) render.polygonMode(MODE_LINE);
        else                render.polygonMode(MODE_FILL);

        cleardevice();
        render.useProgram(shaders[selectShader]);
        render.drawElements();

        outtextxy(WIDTH / 2, 10, shaderName[selectShader]);
        swprintf_s(buf, _T("  vertices   : %zu"), vertices);
        outtextxy(0, 0, buf);
        swprintf_s(buf, _T("  triangles  : %zu"), triangles);
        outtextxy(0, 16, buf);
        outtextxy(0, 32, L"  Scroll       : fov");
        outtextxy(0, 48, L"  Drag LButton    : Rotate");
        outtextxy(0, 64, L"  LButton + Scroll: Forward / Backward");
        outtextxy(0, 80, L"  Press [C]: Color texture");
        outtextxy(0, 96, L"  Press [N]: Normal texture");
        outtextxy(0, 112, L"  Press [B]: Bump texture");
        outtextxy(0, 128, L"  Press [L]: Line Mode");
        outtextxy(0, 144, L"  Press [K]: Cull Face");
        outtextxy(0, 160, L"  Press [Space] : Shaders");
        swprintf_s(buf, _T("  Camera Position: (%.3f, %.3f, %.3f)"), camera.position().x, camera.position().y, camera.position().z);
        outtextxy(0, 176, buf);
        swprintf_s(buf, _T("  Camera Fov: %.3f"), camera.fov());
        outtextxy(0, 192, buf);
        swprintf_s(buf, _T("fps  : %d"), (int)(1.f / deltaTime));
        outtextxy(WIDTH - 96, 0, buf);
        swprintf_s(buf, _T("time: %.1fms"), 1000 * deltaTime);
        outtextxy(WIDTH - 96, 16, buf);
        FlushBatchDraw();
    }
    return 0;
}