//[header]
// This program illustrates how the concept of vector and matrix can be implemented
// in C++. This is a light version of the implementation. It contains the most
// essential methods to manipulate vectors and matrices. It should be enough
// for most projects. Vectors and matrices are really the alphabet as we said
// in the lesson of any graphics application. It's really important you feel
// confortable with these techniques especially with the concepts of
// normalizing vectors, computing their length, computing the dot and cross products
// of two vectors, and the point- and vector-matrix multiplication (and knowing
// the difference between the two).
//[/header]
//[compile]
// c++ geometry.cpp -o geometry -std=c++11
//[/compile]
//[ignore]
// Copyright (C) 2012 www.scratchapixel.com
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//[/ignore]
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <iomanip>
//[comment]
// Implementation of a generic vector class - it will be used to deal with 3D points, vectors and normals.
// The class is implemented as a template. While it may complicate the code a bit, it gives us
// the flexibility later, to specialize the type of the coordinates into anything we want.
// For example: Vec3f if we want the coordinates to be floats or Vec3i if we want the coordinates to be integers.
//
// Vec3 is a standard/common way of naming vectors, points, etc. The OpenEXR and Autodesk libraries
// use this convention for instance.
//[/comment]
template<typename T>
class Vec3
{
public:
Vec3() : x(0), y(0), z(0) {}
Vec3(T xx) : x(xx), y(xx), z(xx) {}
Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {}
Vec3 operator + (const Vec3 &v) const
{ return Vec3(x + v.x, y + v.y, z + v.z); }
Vec3 operator - (const Vec3 &v) const
{ return Vec3(x - v.x, y - v.y, z - v.z); }
Vec3 operator * (const T &r) const
{ return Vec3(x * r, y * r, z * r); }
T dotProduct(const Vec3<T> &v) const
{ return x * v.x + y * v.y + z * v.z; }
T crossProduct(const Vec3<T> &v) const
{ return Vec3<T>(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x); }
T norm() const
{ return x * x + y * y + z * z; }
T length() const
{ return sqrt(norm()); }
//[comment]
// The next two operators are sometimes called access operators or
// accessors. The Vec coordinates can be accessed that way v[0], v[1], v[2],
// rather than using the more traditional form v.x, v.y, v.z. This useful
// when vectors are used in loops: the coordinates can be accessed with the
// loop index (e.g. v[i]).
//[/comment]
const T& operator [] (uint8_t i) const { return (&x)[i]; }
T& operator [] (uint8_t i) { return (&x)[i]; }
Vec3& normalize()
{
T n = norm();
if (n > 0) {
T factor = 1 / sqrt(n);
x *= factor, y *= factor, z *= factor;
}
return *this;
}
friend std::ostream& operator << (std::ostream &s, const Vec3<T> &v)
{
return s << '(' << v.x << ' ' << v.y << ' ' << v.z << ')';
}
T x, y, z;
};
//[comment]
// Now you can specialize the class. We are just showing two examples here. In your code
// you can declare a vector either that way: Vec3<float> a, or that way: Vec3f a
//[/comment]
typedef Vec3<float> Vec3f;
typedef Vec3<int> Vec3i;
//[comment]
// Implementation of a generic 4x4 Matrix class - Same thing here than with the Vec3 class. It uses
// a template which is maybe less useful than with vectors but it can be used to
// define the coefficients of the matrix to be either floats (the most case) or doubles depending
// on our needs.
//
// To use you can either write: Matrix44<float> m; or: Matrix44f m;
//[/comment]
template<typename T>
class Matrix44
{
public:
T x[4][4] = {{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
Matrix44() {}
Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
T i, T j, T k, T l, T m, T n, T o, T p)
{
x[0][0] = a;
x[0][1] = b;
x[0][2] = c;
x[0][3] = d;
x[1][0] = e;
x[1][1] = f;
x[1][2] = g;
x[1][3] = h;
x[2][0] = i;
x[2][1] = j;
x[2][2] = k;
x[2][3] = l;
x[3][0] = m;
x[3][1] = n;
x[3][2] = o;
x[3][3] = p;
}
const T* operator [] (uint8_t i) const { return x[i]; }
T* operator [] (uint8_t i) { return x[i]; }
// Multiply the current matrix with another matrix (rhs)
Matrix44 operator * (const Matrix44& v) const
{
Matrix44 tmp;
multiply (*this, v, tmp);
return tmp;
}
//[comment]
// To make it easier to understand how a matrix multiplication works, the fragment of code
// included within the #if-#else statement, show how this works if you were to iterate
// over the coefficients of the resulting matrix (a). However you will often see this
// multiplication being done using the code contained within the #else-#end statement.
// It is exactly the same as the first fragment only we have litteraly written down
// as a series of operations what would actually result from executing the two for() loops
// contained in the first fragment. It is supposed to be faster, however considering
// matrix multiplicatin is not necessarily that common, this is probably not super
// useful nor really necessary (but nice to have -- and it gives you an example of how
// it can be done, as this how you will this operation implemented in most libraries).
//[/comment]
static void multiply(const Matrix44<T> &a, const Matrix44& b, Matrix44 &c)
{
#if 0
for (uint8_t i = 0; i < 4; ++i) {
for (uint8_t j = 0; j < 4; ++j) {
c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] +
a[i][2] * b[2][j] + a[i][3] * b[3][j];
}
}
#else
// A restric qualified pointer (or reference) is basically a promise
// to the compiler that for the scope of the pointer, the target of the
// pointer will only be accessed through that pointer (and pointers
// copied from it.
const T * __restrict ap = &a.x[0][0];
const T * __restrict bp = &b.x[0][0];
T * __restrict cp = &c.x[0][0];
T a0, a1, a2, a3;
a0 = ap[0];
a1 = ap[1];
a2 = ap[2];
a3 = ap[3];
cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
a0 = ap[4];
a1 = ap[5];
a2 = ap[6];
a3 = ap[7];
cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
a0 = ap[8];
a1 = ap[9];
a2 = ap[10];
a3 = ap[11];
cp[8] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
cp[9] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
a0 = ap[12];
a1 = ap[13];
a2 = ap[14];
a3 = ap[15];
cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
#endif
}
// \brief return a transposed copy of the current matrix as a new matrix
Matrix44 transposed() const
{
#if 0
Matrix44 t;
for (uint8_t i = 0; i < 4; ++i) {
for (uint8_t j = 0; j < 4; ++j) {
t[i][j] = x[j][i];
}
}
return t;
#else
return Matrix44 (x[0][0],
x[1][0],
x[2][0],
x[3][0],
x[0][1],
x[1][1],
x[2][1],
x[3][1],
x[0][2],
x[1][2],
x[2][2],
x[3][2],
x[0][3],
x[1][3],
x[2][3],
x[3][3]);
#endif
}
// \brief transpose itself
Matrix44& transpose ()
{
Matrix44 tmp (x[0][0],
x[1][0],
x[2][0],
x[3][0],
x[0][1],
x[1][1],
x[2][1],
x[3][1],
x[0][2],
x[1][2],
x[2][2],
x[3][2],
x[0][3],
x[1][3],
x[2][3],
x[3][3]);
*this = tmp;
return *this;
}
//[comment]
// This method needs to be used for point-matrix multiplication. Keep in mind
// we don't make the distinction between points and vectors at least from
// a programming point of view, as both (as well as normals) are declared as Vec3.
// However, mathematically they need to be treated differently. Points can be translated
// when translation for vectors is meaningless. Furthermore, points are implicitly
// be considered as having homogeneous coordinates. Thus the w coordinates needs
// to be computed and to convert the coordinates from homogeneous back to Cartesian
// coordinates, we need to divided x, y z by w.
//
// The coordinate w is more often than not equals to 1, but it can be different than
// 1 especially when the matrix is projective matrix (perspective projection matrix).
//[/comment]
template<typename S>
void multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const
{
S a, b, c, w;
a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0];
b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1];
c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2];
w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3];
dst.x = a / w;
dst.y = b / w;
dst.z = c / w;
}
//[comment]
// This method needs to be used for vector-matrix multiplication. Look at the differences
// with the previous method (to compute a point-matrix multiplication). We don't use
// the coefficients in the matrix that account for translation (x[3][0], x[3][1], x[3][2])
// and we don't compute w.
//[/comment]
template<typename S>
void multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const
{
S a, b, c;
a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0];
b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1];
c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2];
dst.x = a;
dst.y = b;
dst.z = c;
}
//[comment]
// Compute the inverse of the matrix using the Gauss-Jordan (or reduced row) elimination method.
// We didn't explain in the lesson on Geometry how the inverse of matrix can be found. Don't
// worry at this point if you don't understand how this works. But we will need to be able to
// compute the inverse of matrices in the first lessons of the "Foundation of 3D Rendering" section,
// which is why we've added this code. For now, you can just use it and rely on it
// for doing what it's supposed to do. If you want to learn how this works though, check the lesson
// on called Matrix Inverse in the "Mathematics and Physics of Computer Graphics" section.
//[/comment]
Matrix44 inverse()
{
int i, j, k;
Matrix44 s;
Matrix44 t (*this);
// Forward elimination
for (i = 0; i < 3 ; i++) {
int pivot = i;
T pivotsize = t[i][i];
if (pivotsize < 0)
pivotsize = -pivotsize;
for (j = i + 1; j < 4; j++) {
T tmp = t[j][i];
if (tmp < 0)
tmp = -tmp;
if (tmp > pivotsize) {
pivot = j;
pivotsize = tmp;
}
}
if (pivotsize == 0) {
// Cannot invert singular matrix
return Matrix44();
}
if (pivot != i) {
for (j = 0; j < 4; j++) {
T tmp;
tmp = t[i][j];
t[i][j] = t[pivot][j];
t[pivot][j] = tmp;
tmp = s[i][j];
s[i][j] = s[pivot][j];
s[pivot][j] = tmp;
}
}
for (j = i + 1; j < 4; j++) {
T f = t[j][i] / t[i][i];
for (k = 0; k < 4; k++) {
t[j][k] -= f * t[i][k];
s[j][k] -= f * s[i][k];
}
}
}
// Backward substitution
for (i = 3; i >= 0; --i) {
T f;
if ((f = t[i][i]) == 0) {
// Cannot invert singular matrix
return Matrix44();
}
for (j = 0; j < 4; j++) {
t[i][j] /= f;
s[i][j] /= f;
}
for (j = 0; j < i; j++) {
f = t[j][i];
for (k = 0; k < 4; k++) {
t[j][k] -= f * t[i][k];
s[j][k] -= f * s[i][k];
}
}
}
return s;
}
// \brief set current matrix to its inverse
const Matrix44<T>& invert()
{
*this = inverse();
return *this;
}
friend std::ostream& operator << (std::ostream &s, const Matrix44 &m)
{
std::ios_base::fmtflags oldFlags = s.flags();
int width = 12; // total with of the displayed number
s.precision(5); // control the number of displayed decimals
s.setf (std::ios_base::fixed);
s << "(" << std::setw (width) << m[0][0] <<
" " << std::setw (width) << m[0][1] <<
" " << std::setw (width) << m[0][2] <<
" " << std::setw (width) << m[0][3] << "\n" <<
" " << std::setw (width) << m[1][0] <<
" " << std::setw (width) << m[1][1] <<
" " << std::setw (width) << m[1][2] <<
" " << std::setw (width) << m[1][3] << "\n" <<
" " << std::setw (width) << m[2][0] <<
" " << std::setw (width) << m[2][1] <<
" " << std::setw (width) << m[2][2] <<
" " << std::setw (width) << m[2][3] << "\n" <<
" " << std::setw (width) << m[3][0] <<
" " << std::setw (width) << m[3][1] <<
" " << std::setw (width) << m[3][2] <<
" " << std::setw (width) << m[3][3] << ")\n";
s.flags (oldFlags);
return s;
}
};
typedef Matrix44<float> Matrix44f;
//[comment]
// Testing our code. To test the matrix inversion code, we used Maya to output
// the values of a matrix and its inverse (check the video at the top of this page). Of course this implies
// that Maya actually does the right thing, but we can probably agree, that is actually does;).
// These are the values for the input matrix:
//
// 0.707107 0 -0.707107 0 -0.331295 0.883452 -0.331295 0 0.624695 0.468521 0.624695 0 4.000574 3.00043 4.000574 1
//
// Given the input matrix, the inverse matrix computed by our code should match the following values:
//
// 0.707107 -0.331295 0.624695 0 0 0.883452 0.468521 0 -0.707107 -0.331295 0.624695 0 0 0 -6.404043 1
//[/comment]
int main(int argc, char **argv)
{
Vec3f v(0, 1, 2);
std::cerr << v << std::endl;
Matrix44f a, b, c;
c = a * b;
Matrix44f d(0.707107, 0, -0.707107, 0, -0.331295, 0.883452, -0.331295, 0, 0.624695, 0.468521, 0.624695, 0, 4.000574, 3.00043, 4.000574, 1);
std::cerr << d << std::endl;
d.invert();
std::cerr << d << std::endl;
return 0;
}
零基础学图形学(14) 几何知识——Geometry源码
最新推荐文章于 2022-11-09 23:09:37 发布