I believe I must be able to beat 130 FPS.
...even Quake2 was running faster
I must admit that I got my inspiration from OpenGL – resembling a few parts of it (at least, how I believe it could work) and with some simplifications.
I must admit also that I didn't understand what exactly the caculateBarycentricWeights() function is good for.
However, I based my implementation on the following concept:
rasterizer must work in screen space
3D vertices have to be converted to screen space before rasterizing
The inner rastering loop has to be as simple as possible – branches (i.e. if and a complex function) in the most inner loop are counter-productive for sure.
Hence, I made a new rasterize() function. Thereby, I considered that rasterize() has to fill horizontal screen lines. For this, I sort the 3 triangle vertices by their y components. (Note, that vertices have to be transformed to screen space beforehand.) In regular case, this allows to cut the triangle horizontally into two parts:
- the upper with peak above of horizontal base
- the lower with peak below of horizontal base.
The source code follows below but the sketch may help to make all the interpolations a little bit more understandable:
![sketch of rasterize()]()
To convert the vtcs beforehand to screen space, I used two matrices:
const Mat4x4f matProj(InitScale, 1.0f / Width, 1.0f / Height, 1.0f);
const Mat4x4f matScreen
= Mat4x4f(InitScale, 0.5f * Width, 0.5f * Height, 1.0f)
* Mat4x4f(InitTrans, Vec3f(1.0f, 1.0f, 0.0f));
const Mat4x4f mat = matScreen * matProj /* * matView * matModel */;
where
matProj is the projection matrix to scale model coordinates to clip spacematScreen is responsible to convert clip space to screen space.
Clip space is a space where visible part is in range [-1, 1] for x, y, and z direction.
Screen space has ranges [0, width) and [0, height) where I used width = 1024 and height = 768 according to OPs requirement.
When I create OpenGL programs it's (some kind of) annoying tradition that I everytimes start with a blue screen (as blue is my favorite clear color). This is mostly caused by confusing some of the matrix operations and followed by the effort to find and fix these operations. Therefore, I tried to visualize the rasterization result to become sure that my benchmarking is not overly enthusiastic due to the fact that simply any part of triangle is outside of view and hence clipped away:
![snapshot of demo using the transformations of benchmark]()
I would count this a visual proof. To achieve this, I wrapped the rasterizer sample into a simple Qt application where the FBO data is stored in a QImage which in turn is applied to a QPixmap displayed with a QLabel. Once I got this, I thought it would be nice to add some kind of simple animation by modifying the model-view matrix over time. So, I got the following sample:
#include <cstdint>
#include <algorithm>
#include <chrono>
#include <iostream>
#include <vector>
#include "linmath.h"
typedef unsigned uint;
typedef std::uint32_t uint32;
const int Width = 1024, Height = 768;
class FBO {
public:
const int width, height;
private:
std::vector<uint32> _rgba; // pixel buffer
public:
FBO(int width, int height):
width(width), height(height),
_rgba(width * height, 0)
{ }
~FBO() = default;
FBO(const FBO&) = delete;
FBO& operator=(const FBO&) = delete;
void clear(uint32 rgba) { std::fill(_rgba.begin(), _rgba.end(), rgba); }
void set(int x, int y, uint32 rgba)
{
const size_t i = y * width + x;
_rgba[i] = rgba;
}
size_t getI(int y) const { return y * width; }
size_t getI(int x, int y) const { return y * width + x; }
void set(size_t i, uint32 rgba) { _rgba[i] = rgba; }
const std::vector<uint32>& getRGBA() const { return _rgba; }
};
void rasterize(FBO &fbo, const Vec3f vtcs[3], const uint32 rgba)
{
// sort vertices by y coordinates
uint iVtcs[3] = { 0, 1, 2 };
if (vtcs[iVtcs[0]].y > vtcs[iVtcs[1]].y) std::swap(iVtcs[0], iVtcs[1]);
if (vtcs[iVtcs[1]].y > vtcs[iVtcs[2]].y) std::swap(iVtcs[1], iVtcs[2]);
if (vtcs[iVtcs[0]].y > vtcs[iVtcs[1]].y) std::swap(iVtcs[0], iVtcs[1]);
const Vec3f vtcsS[3] = { vtcs[iVtcs[0]], vtcs[iVtcs[1]], vtcs[iVtcs[2]] };
const float yM = vtcs[1].y;
// cut triangle in upper and lower part
const float xT = vtcsS[0].x;
float xML = vtcsS[1].x;
const float f = (yM - vtcsS[0].y) / (vtcsS[2].y - vtcsS[0].y);
float xMR = (1.0f - f) * xT + f * vtcsS[2].x;
if (xML > xMR) std::swap(xML, xMR);
// draw upper part of triangle
if (vtcs[iVtcs[0]].y < yM) {
const float dY = yM - vtcsS[0].y;
for (int y = std::max((int)vtcsS[0].y, 0),
yE = std::min((int)(yM + 0.5f), fbo.height);
y < yE; ++y) {
const float f1 = (yM - y) / dY, f0 = 1.0f - f1;
const float xL = f0 * xT + f1 * xML, xR = f0 * xT + f1 * xMR;
const size_t i = fbo.getI(y);
for (int x = std::max((int)xL, 0),
xE = std::min((int)(xR + 0.5f), fbo.width);
x < xE; ++x) {
fbo.set(i + x, rgba);
}
}
} // else upper edge horizontal
// draw lower part of triangle
if (yM < vtcs[2].y) {
const float xB = vtcsS[2].x;
const float dY = vtcsS[2].y - yM;
for (int y = std::max((int)yM, 0),
yE = std::min((int)(vtcsS[2].y + 0.5f), fbo.height);
y < yE; ++y) {
const float f1 = (y - yM) / dY, f0 = 1.0f - f1;
const float xL = f0 * xML + f1 * xB, xR = f0 * xMR + f1 * xB;
const size_t i = fbo.getI(y);
for (int x = std::max((int)xL, 0),
xE = std::min((int)(xR + 0.5f), fbo.width);
x < xE; ++x) {
fbo.set(i + x, rgba);
}
}
} // else lower edge horizontal
}
template <typename VALUE>
Vec3T<VALUE> transformPoint(const Mat4x4T<VALUE> &mat, const Vec3T<VALUE> &pt)
{
Vec4T<VALUE> pt_ = mat * Vec4T<VALUE>(pt.x, pt.y, pt.z, (VALUE)1);
return pt_.w != (VALUE)0
? Vec3T<VALUE>(pt_.x / pt_.w, pt_.y / pt_.w, pt_.z / pt_.w)
: Vec3T<VALUE>(pt_.x, pt_.y, pt_.z);
}
typedef std::chrono::high_resolution_clock HiResClock;
typedef std::chrono::microseconds MicroSecs;
int mainBench()
{
const Mat4x4f matProj(InitScale, 1.0f / Width, 1.0f / Height, 1.0f);
const Mat4x4f matScreen
= Mat4x4f(InitScale, 0.5f * Width, 0.5f * Height, 1.0f)
* Mat4x4f(InitTrans, Vec3f(1.0f, 1.0f, 0.0f));
const Mat4x4f mat = matScreen * matProj /* * matView * matModel */;
FBO fbo(Width, Height);
HiResClock::time_point start = HiResClock::now();
size_t fps = 0;
for (;;) {
static Vec3f v0 = { -2000.0f, -2000.0, 10.0 };
static Vec3f v1 = { 2000.0f, -2000.0, 10.0 };
static Vec3f v2 = { 0.0f, 2000.0, 10.0 };
const Vec3f vtcs[] = {
transformPoint(mat, v0), transformPoint(mat, v1), transformPoint(mat, v2)
};
rasterize(fbo, vtcs, 0xff0000ff);
++fps;
HiResClock::time_point now = HiResClock::now();
auto timeDiff
= std::chrono::duration_cast<MicroSecs>(now - start).count();
static const long oneSecond = 1000000;
if (timeDiff >= oneSecond) {
std::cout << fps << " " << std::flush;
fps = 0;
start = now;
}
}
}
#include <stack>
#include <QtWidgets>
struct RenderContext {
FBO fbo; // frame buffer object
Mat4x4f matModelView;
Mat4x4f matProj;
RenderContext(int width, int height):
fbo(width, height),
matModelView(InitIdent),
matProj(InitIdent)
{ }
~RenderContext() = default;
};
void render(RenderContext &context)
{
context.fbo.clear(0xffffcc88);
static Vec3f v0 = { -2000.0f, -2000.0, 10.0 };
static Vec3f v1 = { 2000.0f, -2000.0, 10.0 };
static Vec3f v2 = { 0.0f, 2000.0, 10.0 };
#if 0 // test transformations of mainBench()
const Mat4x4f matProj(InitScale, 1.0f / Width, 1.0f / Height, 1.0f);
const Mat4x4f matScreen
= Mat4x4f(InitScale, 0.5f * Width, 0.5f * Height, 1.0f)
* Mat4x4f(InitTrans, Vec3f(1.0f, 1.0f, 0.0f));
const Mat4x4f mat = matScreen * matProj /* * matView * matModel */;
#else // make a simple animation
const Mat4x4f matScreen
= Mat4x4f(InitScale, 0.5f * context.fbo.width, 0.5f * context.fbo.height, 1.0f)
* Mat4x4f(InitTrans, Vec3f(1.0f, 1.0f, 0.0f));
const Mat4x4f mat = matScreen * context.matProj * context.matModelView;
#endif // 0
const Vec3f vtcs[] = {
transformPoint(mat, v0), transformPoint(mat, v1), transformPoint(mat, v2)
};
rasterize(context.fbo, vtcs, 0xff0000ff);
}
int main(int argc, char **argv)
{
// evaluate arguments
const bool gui = argc > 1
&& (strcmp(argv[1], "gui") == 0 || strcmp(argv[1], "-gui") == 0);
// without arguments: start benchmark
if (!gui) return mainBench();
// remove argument "gui"
for (char **arg = argv + 1; *arg; ++arg) arg [0] = arg[1];
// Qt Demo
qDebug() << "Qt Version:" << QT_VERSION_STR;
QApplication app(argc, argv);
RenderContext context(Width, Height);
// setup GUI
QPixmap qPixmapImg(Width, Height);
QLabel qLblImg;
qLblImg.setWindowTitle("Software Rasterizer Demo");
qLblImg.setPixmap(qPixmapImg);
qLblImg.show();
// Qt timer for periodic rendering/animation
QTime qTime(0, 0);
QTimer qTimer;
qTimer.setInterval(50);
// install signal handlers
QObject::connect(&qTimer, &QTimer::timeout,
[&]() {
// simple animation
const float r = 1.0f, t = 0.001f * qTime.elapsed();
context.matModelView
= Mat4x4f(InitTrans, Vec3f(r * sin(t), r * cos(t), 0.0f))
* Mat4x4f(InitScale, 0.0001f, 0.0001f, 0.0001f);
// rendering
render(context);
// transfer FBO to qLblImg
const QImage qImg((uchar*)context.fbo.getRGBA().data(),
Width, Height, QImage::Format_RGBA8888);
qPixmapImg.convertFromImage(qImg);
qLblImg.setPixmap(qPixmapImg);
});
// runtime loop
qTime.start(); qTimer.start();
return app.exec();
}
Instead of gmtl (which I don't know), I used linmath.h (left from earlier MCVEs I wrote in the past):
#ifndef LIN_MATH_H
#define LIN_MATH_H
#include <iostream>
#include <cassert>
#include <cmath>
double Pi = 3.1415926535897932384626433832795;
template <typename VALUE>
inline VALUE degToRad(VALUE angle)
{
return (VALUE)Pi * angle / (VALUE)180;
}
template <typename VALUE>
inline VALUE radToDeg(VALUE angle)
{
return (VALUE)180 * angle / (VALUE)Pi;
}
template <typename VALUE>
struct Vec2T {
VALUE x, y;
Vec2T() { }
Vec2T(VALUE x, VALUE y): x(x), y(y) { }
};
template <typename VALUE>
VALUE length(const Vec2T<VALUE> &vec)
{
return sqrt(vec.x * vec.x + vec.y * vec.y);
}
template <typename VALUE>
std::ostream& operator<<(std::ostream &out, const Vec2T<VALUE> &v)
{
return out << "( " << v.x << ", " << v.y << " )";
}
typedef Vec2T<float> Vec2f;
typedef Vec2T<double> Vec2;
template <typename VALUE>
struct Vec3T {
VALUE x, y, z;
Vec3T() { }
Vec3T(VALUE x, VALUE y, VALUE z): x(x), y(y), z(z) { }
Vec3T(const Vec2T<VALUE> &xy, VALUE z): x(xy.x), y(xy.y), z(z) { }
explicit operator Vec2T<VALUE>() const { return Vec2T<VALUE>(x, y); }
};
template <typename VALUE>
std::ostream& operator<<(std::ostream &out, const Vec3T<VALUE> &v)
{
return out << "( " << v.x << ", " << v.y << ", " << v.z << " )";
}
typedef Vec3T<float> Vec3f;
typedef Vec3T<double> Vec3;
template <typename VALUE>
struct Vec4T {
VALUE x, y, z, w;
Vec4T() { }
Vec4T(VALUE x, VALUE y, VALUE z, VALUE w): x(x), y(y), z(z), w(w) { }
Vec4T(const Vec2T<VALUE> &xy, VALUE z, VALUE w):
x(xy.x), y(xy.y), z(z), w(w)
{ }
Vec4T(const Vec3T<VALUE> &xyz, VALUE w):
x(xyz.x), y(xyz.y), z(xyz.z), w(w)
{ }
explicit operator Vec2T<VALUE>() const { return Vec2T<VALUE>(x, y); }
explicit operator Vec3T<VALUE>() const { return Vec3T<VALUE>(x, y, z); }
};
template <typename VALUE>
std::ostream& operator<<(std::ostream &out, const Vec4T<VALUE> &v)
{
return out << "( " << v.x << ", " << v.y << ", " << v.z << ", " << v.w << " )";
}
typedef Vec4T<float> Vec4f;
typedef Vec4T<double> Vec4;
enum ArgInitIdent { InitIdent };
enum ArgInitTrans { InitTrans };
enum ArgInitRot { InitRot };
enum ArgInitRotX { InitRotX };
enum ArgInitRotY { InitRotY };
enum ArgInitRotZ { InitRotZ };
enum ArgInitScale { InitScale };
template <typename VALUE>
struct Mat4x4T {
union {
VALUE comp[4 * 4];
struct {
VALUE _00, _01, _02, _03;
VALUE _10, _11, _12, _13;
VALUE _20, _21, _22, _23;
VALUE _30, _31, _32, _33;
};
};
// constructor to build a matrix by elements
Mat4x4T(
VALUE _00, VALUE _01, VALUE _02, VALUE _03,
VALUE _10, VALUE _11, VALUE _12, VALUE _13,
VALUE _20, VALUE _21, VALUE _22, VALUE _23,
VALUE _30, VALUE _31, VALUE _32, VALUE _33):
_00(_00), _01(_01), _02(_02), _03(_03),
_10(_10), _11(_11), _12(_12), _13(_13),
_20(_20), _21(_21), _22(_22), _23(_23),
_30(_30), _31(_31), _32(_32), _33(_33)
{ }
// constructor to build an identity matrix
Mat4x4T(ArgInitIdent):
_00((VALUE)1), _01((VALUE)0), _02((VALUE)0), _03((VALUE)0),
_10((VALUE)0), _11((VALUE)1), _12((VALUE)0), _13((VALUE)0),
_20((VALUE)0), _21((VALUE)0), _22((VALUE)1), _23((VALUE)0),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{ }
// constructor to build a matrix for translation
Mat4x4T(ArgInitTrans, const Vec3T<VALUE> &t):
_00((VALUE)1), _01((VALUE)0), _02((VALUE)0), _03((VALUE)t.x),
_10((VALUE)0), _11((VALUE)1), _12((VALUE)0), _13((VALUE)t.y),
_20((VALUE)0), _21((VALUE)0), _22((VALUE)1), _23((VALUE)t.z),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{ }
// constructor to build a matrix for rotation about axis
Mat4x4T(ArgInitRot, const Vec3T<VALUE> &axis, VALUE angle):
_03((VALUE)0), _13((VALUE)0), _23((VALUE)0),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{
//axis.normalize();
const VALUE sinAngle = sin(angle), cosAngle = cos(angle);
const VALUE xx = axis.x * axis.x, xy = axis.x * axis.y;
const VALUE xz = axis.x * axis.z, yy = axis.y * axis.y;
const VALUE yz = axis.y * axis.z, zz = axis.z * axis.z;
_00 = xx + cosAngle * ((VALUE)1 - xx) /* + sinAngle * 0 */;
_01 = xy - cosAngle * xy - sinAngle * axis.z;
_02 = xz - cosAngle * xz + sinAngle * axis.y;
_10 = xy - cosAngle * xy + sinAngle * axis.z;
_11 = yy + cosAngle * ((VALUE)1 - yy) /* + sinAngle * 0 */;
_12 = yz - cosAngle * yz - sinAngle * axis.x;
_20 = xz - cosAngle * xz - sinAngle * axis.y;
_21 = yz - cosAngle * yz + sinAngle * axis.x;
_22 = zz + cosAngle * ((VALUE)1 - zz) /* + sinAngle * 0 */;
}
// constructor to build a matrix for rotation about x axis
Mat4x4T(ArgInitRotX, VALUE angle):
_00((VALUE)1), _01((VALUE)0), _02((VALUE)0), _03((VALUE)0),
_10((VALUE)0), _11(cos(angle)), _12(-sin(angle)), _13((VALUE)0),
_20((VALUE)0), _21(sin(angle)), _22(cos(angle)), _23((VALUE)0),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{ }
// constructor to build a matrix for rotation about y axis
Mat4x4T(ArgInitRotY, VALUE angle):
_00(cos(angle)), _01((VALUE)0), _02(sin(angle)), _03((VALUE)0),
_10((VALUE)0), _11((VALUE)1), _12((VALUE)0), _13((VALUE)0),
_20(-sin(angle)), _21((VALUE)0), _22(cos(angle)), _23((VALUE)0),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{ }
// constructor to build a matrix for rotation about z axis
Mat4x4T(ArgInitRotZ, VALUE angle):
_00(cos(angle)), _01(-sin(angle)), _02((VALUE)0), _03((VALUE)0),
_10(sin(angle)), _11(cos(angle)), _12((VALUE)0), _13((VALUE)0),
_20((VALUE)0), _21((VALUE)0), _22((VALUE)1), _23((VALUE)0),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{ }
// constructor to build a matrix for scaling
Mat4x4T(ArgInitScale, VALUE sx, VALUE sy, VALUE sz):
_00((VALUE)sx), _01((VALUE)0), _02((VALUE)0), _03((VALUE)0),
_10((VALUE)0), _11((VALUE)sy), _12((VALUE)0), _13((VALUE)0),
_20((VALUE)0), _21((VALUE)0), _22((VALUE)sz), _23((VALUE)0),
_30((VALUE)0), _31((VALUE)0), _32((VALUE)0), _33((VALUE)1)
{ }
// operator to allow access with [][]
double* operator [] (int i)
{
assert(i >= 0 && i < 4);
return comp + 4 * i;
}
// operator to allow access with [][]
const double* operator [] (int i) const
{
assert(i >= 0 && i < 4);
return comp + 4 * i;
}
// multiply matrix with matrix -> matrix
Mat4x4T operator * (const Mat4x4T &mat) const
{
return Mat4x4T(
_00 * mat._00 + _01 * mat._10 + _02 * mat._20 + _03 * mat._30,
_00 * mat._01 + _01 * mat._11 + _02 * mat._21 + _03 * mat._31,
_00 * mat._02 + _01 * mat._12 + _02 * mat._22 + _03 * mat._32,
_00 * mat._03 + _01 * mat._13 + _02 * mat._23 + _03 * mat._33,
_10 * mat._00 + _11 * mat._10 + _12 * mat._20 + _13 * mat._30,
_10 * mat._01 + _11 * mat._11 + _12 * mat._21 + _13 * mat._31,
_10 * mat._02 + _11 * mat._12 + _12 * mat._22 + _13 * mat._32,
_10 * mat._03 + _11 * mat._13 + _12 * mat._23 + _13 * mat._33,
_20 * mat._00 + _21 * mat._10 + _22 * mat._20 + _23 * mat._30,
_20 * mat._01 + _21 * mat._11 + _22 * mat._21 + _23 * mat._31,
_20 * mat._02 + _21 * mat._12 + _22 * mat._22 + _23 * mat._32,
_20 * mat._03 + _21 * mat._13 + _22 * mat._23 + _23 * mat._33,
_30 * mat._00 + _31 * mat._10 + _32 * mat._20 + _33 * mat._30,
_30 * mat._01 + _31 * mat._11 + _32 * mat._21 + _33 * mat._31,
_30 * mat._02 + _31 * mat._12 + _32 * mat._22 + _33 * mat._32,
_30 * mat._03 + _31 * mat._13 + _32 * mat._23 + _33 * mat._33);
}
// multiply matrix with vector -> vector
Vec4T<VALUE> operator * (const Vec4T<VALUE> &vec) const
{
return Vec4T<VALUE>(
_00 * vec.x + _01 * vec.y + _02 * vec.z + _03 * vec.w,
_10 * vec.x + _11 * vec.y + _12 * vec.z + _13 * vec.w,
_20 * vec.x + _21 * vec.y + _22 * vec.z + _23 * vec.w,
_30 * vec.x + _31 * vec.y + _32 * vec.z + _33 * vec.w);
}
};
template <typename VALUE>
std::ostream& operator<<(std::ostream &out, const Mat4x4T<VALUE> &m)
{
return out
<< m._00 << '\t' << m._01 << '\t' << m._02 << '\t' << m._03 << '\n'
<< m._10 << '\t' << m._11 << '\t' << m._12 << '\t' << m._13 << '\n'
<< m._20 << '\t' << m._21 << '\t' << m._22 << '\t' << m._23 << '\n'
<< m._30 << '\t' << m._31 << '\t' << m._32 << '\t' << m._33 << '\n';
}
typedef Mat4x4T<float> Mat4x4f;
typedef Mat4x4T<double> Mat4x4;
// enumeration of rotation axes
enum RotAxis {
RotX, // rotation about x axis
RotY, // rotation about y axis
RotZ // rotation about z axis
};
// enumeration of possible Euler angles
enum EulerAngle {
RotXYX = RotX + 3 * RotY + 9 * RotX, // 0 + 3 + 0 = 3
RotXYZ = RotX + 3 * RotY + 9 * RotZ, // 0 + 3 + 18 = 21
RotXZX = RotX + 3 * RotZ + 9 * RotX, // 0 + 6 + 0 = 6
RotXZY = RotX + 3 * RotZ + 9 * RotY, // 0 + 6 + 9 = 15
RotYXY = RotY + 3 * RotX + 9 * RotY, // 1 + 0 + 9 = 10
RotYXZ = RotY + 3 * RotX + 9 * RotZ, // 1 + 0 + 18 = 19
RotYZX = RotY + 3 * RotZ + 9 * RotX, // 1 + 6 + 0 = 7
RotYZY = RotY + 3 * RotZ + 9 * RotY, // 1 + 6 + 9 = 16
RotZXY = RotZ + 3 * RotX + 9 * RotY, // 2 + 0 + 9 = 11
RotZXZ = RotZ + 3 * RotX + 9 * RotZ, // 2 + 0 + 18 = 20
RotZYX = RotZ + 3 * RotY + 9 * RotX, // 2 + 3 + 0 = 5
RotZYZ = RotZ + 3 * RotY + 9 * RotZ, // 2 + 3 + 18 = 23
RotHPR = RotZXY, // used in OpenGL Performer
RotABC = RotZYX // used in German engineering
};
/* decomposes the combined EULER angle type into the corresponding
* individual EULER angle axis types.
*/
inline void decompose(
EulerAngle type, RotAxis &axis1, RotAxis &axis2, RotAxis &axis3)
{
unsigned type_ = (unsigned)type;
axis1 = (RotAxis)(type_ % 3); type_ /= 3;
axis2 = (RotAxis)(type_ % 3); type_ /= 3;
axis3 = (RotAxis)type_;
}
template <typename VALUE>
Mat4x4T<VALUE> makeEuler(
EulerAngle mode, VALUE rot1, VALUE rot2, VALUE rot3)
{
RotAxis axis1, axis2, axis3;
decompose(mode, axis1, axis2, axis3);
const static VALUE axes[3][3] = {
{ (VALUE)1, (VALUE)0, (VALUE)0 },
{ (VALUE)0, (VALUE)1, (VALUE)0 },
{ (VALUE)0, (VALUE)0, (VALUE)1 }
};
return
Mat4x4T<VALUE>(InitRot,
Vec3T<VALUE>(axes[axis1][0], axes[axis1][1], axes[axis1][2]),
rot1)
* Mat4x4T<VALUE>(InitRot,
Vec3T<VALUE>(axes[axis2][0], axes[axis2][1], axes[axis2][2]),
rot2)
* Mat4x4T<VALUE>(InitRot,
Vec3T<VALUE>(axes[axis3][0], axes[axis3][1], axes[axis3][2]),
rot3);
}
/* decomposes a rotation matrix into EULER angles.
*
* It is necessary that the upper left 3x3 matrix is composed of rotations
* only. Translational parts are not considered.
* Other transformations (e.g. scaling, shearing, projection) may cause
* wrong results.
*/
template <typename VALUE>
void decompose(
const Mat4x4T<VALUE> &mat,
RotAxis axis1, RotAxis axis2, RotAxis axis3,
VALUE &angle1, VALUE &angle2, VALUE &angle3)
{
assert(axis1 != axis2 && axis2 != axis3);
/* This is ported from EulerAngles.h of the Eigen library. */
const int odd = (axis1 + 1) % 3 == axis2 ? 0 : 1;
const int i = axis1;
const int j = (axis1 + 1 + odd) % 3;
const int k = (axis1 + 2 - odd) % 3;
if (axis1 == axis3) {
angle1 = atan2(mat[j][i], mat[k][i]);
if ((odd && angle1 < (VALUE)0) || (!odd && angle1 > (VALUE)0)) {
angle1 = angle1 > (VALUE)0 ? angle1 - (VALUE)Pi : angle1 + (VALUE)Pi;
const VALUE s2 = length(Vec2T<VALUE>(mat[j][i], mat[k][i]));
angle2 = -atan2(s2, mat[i][i]);
} else {
const VALUE s2 = length(Vec2T<VALUE>(mat[j][i], mat[k][i]));
angle2 = atan2(s2, mat[i][i]);
}
const VALUE s1 = sin(angle1);
const VALUE c1 = cos(angle1);
angle3 = atan2(c1 * mat[j][k] - s1 * mat[k][k],
c1 * mat[j][j] - s1 * mat[k][j]);
} else {
angle1 = atan2(mat[j][k], mat[k][k]);
const VALUE c2 = length(Vec2T<VALUE>(mat[i][i], mat[i][j]));
if ((odd && angle1<(VALUE)0) || (!odd && angle1 > (VALUE)0)) {
angle1 = (angle1 > (VALUE)0)
? angle1 - (VALUE)Pi : angle1 + (VALUE)Pi;
angle2 = atan2(-mat[i][k], -c2);
} else angle2 = atan2(-mat[i][k], c2);
const VALUE s1 = sin(angle1);
const VALUE c1 = cos(angle1);
angle3 = atan2(s1 * mat[k][i] - c1 * mat[j][i],
c1 * mat[j][j] - s1 * mat[k][j]);
}
if (!odd) {
angle1 = -angle1; angle2 = -angle2; angle3 = -angle3;
}
}
#endif // LIN_MATH_H
Before I added the Qt code I compiled simply with:
$ g++ -std=c++11 -o rasterizerSimple rasterizerSimple.cc
To compile with Qt, I wrote a Qt project rasterizerSimple.pro:
SOURCES = rasterizerSimple.cc
QT += widgets
Compiled and tested on cygwin64 on Windows 10:
$ qmake-qt5 rasterizerSimple.pro
$ make
$ ./rasterizerSimple
2724 2921 3015 2781 2800 2805 2859 2783 2931 2871 2902 2983 2995 2882 2940 2878 3007 2993 3066 3067 3118 3144 2849 3084 3020 3005 3030 2991 2999 3065 2941 3123 3119 2905 3135 2938
CtrlC
My laptop has an Intel i7 (and produced a bit cooler noise when I let run the program a bit longer).
That looks not that bad – an average of 3000 triangles per second. (It beats the the 130 FPS of OP if his i7 is not much slower than mine.)
May be, the FPS could be pushed a bit more by further improvements. For example, interpolation could be done using Bresenham's line algorithm though it does not appear in most inner loop. To optimize most inner loop, the filling of horizontal line should be done in a faster fashion (although I'm not quite sure how to achieve this).
To start the Qt visualization, the application has to be started with gui:
$ ./rasterizerSimple gui
![snapshot of rasterizerSimple with simple animation and Qt display]()
This sample inspired me to a toy project No-GL 3D Renderer which can be found on github.
