diff --git a/geometry.cpp b/geometry.cpp
index 5d769e28..4324f634 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -10,3 +10,103 @@ template <> template <> Vec3<int>::Vec3<>(const Vec3<float> &v) : x(int(v.x+.5))
 template <> template <> Vec3<float>::Vec3<>(const Vec3<int> &v) : x(v.x), y(v.y), z(v.z) {
 }
 
+
+Matrix::Matrix(int r, int c) : m(std::vector<std::vector<float> >(r, std::vector<float>(c, 0.f))), rows(r), cols(c) { }
+
+int Matrix::nrows() {
+    return rows;
+}
+
+int Matrix::ncols() {
+    return cols;
+}
+
+Matrix Matrix::identity(int dimensions) {
+    Matrix E(dimensions, dimensions);
+    for (int i=0; i<dimensions; i++) {
+        for (int j=0; j<dimensions; j++) {
+            E[i][j] = (i==j ? 1.f : 0.f);
+        }
+    }
+    return E;
+}
+
+std::vector<float>& Matrix::operator[](const int i) {
+    assert(i>=0 && i<rows);
+    return m[i];
+}
+
+Matrix Matrix::operator*(const Matrix& a) {
+    assert(cols == a.rows);
+    Matrix result(rows, a.cols);
+    for (int i=0; i<rows; i++) {
+        for (int j=0; j<a.cols; j++) {
+            result.m[i][j] = 0.f;
+            for (int k=0; k<cols; k++) {
+                result.m[i][j] += m[i][k]*a.m[k][j];
+            }
+        }
+    }
+    return result;
+}
+
+Matrix Matrix::transpose() {
+    Matrix result(cols, rows);
+    for(int i=0; i<rows; i++)
+        for(int j=0; j<cols; j++)
+            result[j][i] = m[i][j];
+    return result;
+}
+
+Matrix Matrix::inverse() {
+    assert(rows==cols);
+    // augmenting the square matrix with the identity matrix of the same dimensions a => [ai]
+    Matrix result(rows, cols*2);
+    for(int i=0; i<rows; i++)
+        for(int j=0; j<cols; j++)
+            result[i][j] = m[i][j];
+    for(int i=0; i<rows; i++)
+        result[i][i+cols] = 1;
+    // first pass
+    for (int i=0; i<rows-1; i++) {
+        // normalize the first row
+        for(int j=result.cols-1; j>=0; j--)
+            result[i][j] /= result[i][i];
+        for (int k=i+1; k<rows; k++) {
+            float coeff = result[k][i];
+            for (int j=0; j<result.cols; j++) {
+                result[k][j] -= result[i][j]*coeff;
+            }
+        }
+    }
+    // normalize the last row
+    for(int j=result.cols-1; j>=rows-1; j--)
+        result[rows-1][j] /= result[rows-1][rows-1];
+    // second pass
+    for (int i=rows-1; i>0; i--) {
+        for (int k=i-1; k>=0; k--) {
+            float coeff = result[k][i];
+            for (int j=0; j<result.cols; j++) {
+                result[k][j] -= result[i][j]*coeff;
+            }
+        }
+    }
+    // cut the identity matrix back
+    Matrix truncate(rows, cols);
+    for(int i=0; i<rows; i++)
+        for(int j=0; j<cols; j++)
+            truncate[i][j] = result[i][j+cols];
+    return truncate;
+}
+
+std::ostream& operator<<(std::ostream& s, Matrix& m) {
+    for (int i=0; i<m.nrows(); i++)  {
+        for (int j=0; j<m.ncols(); j++) {
+            s << m[i][j];
+            if (j<m.ncols()-1) s << "\t";
+        }
+        s << "\n";
+    }
+    return s;
+}
+
diff --git a/geometry.h b/geometry.h
index b5081d34..141d5367 100644
--- a/geometry.h
+++ b/geometry.h
@@ -66,4 +66,28 @@ template <class t> std::ostream& operator<<(std::ostream& s, Vec3<t>& v) {
     return s;
 }
 
+//////////////////////////////////////////////////////////////////////////////////////////////
+
+const int DEFAULT_ALLOC=4;
+
+class Matrix {
+    std::vector<std::vector<float> > m;
+    int rows, cols;
+public:
+    Matrix(int r=DEFAULT_ALLOC, int c=DEFAULT_ALLOC);
+    inline int nrows();
+    inline int ncols();
+
+    static Matrix identity(int dimensions);
+    std::vector<float>& operator[](const int i);
+    Matrix operator*(const Matrix& a);
+    Matrix transpose();
+    Matrix inverse();
+
+    friend std::ostream& operator<<(std::ostream& s, Matrix& m);
+};
+
+/////////////////////////////////////////////////////////////////////////////////////////////
+
+
 #endif //__GEOMETRY_H__
diff --git a/main.cpp b/main.cpp
index db5bd63c..f80d94f9 100644
--- a/main.cpp
+++ b/main.cpp
@@ -12,6 +12,32 @@ const int depth  = 255;
 Model *model = NULL;
 int *zbuffer = NULL;
 Vec3f light_dir(0,0,-1);
+Vec3f camera(0,0,3);
+
+Vec3f m2v(Matrix m) {
+    return Vec3f(m[0][0]/m[3][0], m[1][0]/m[3][0], m[2][0]/m[3][0]);
+}
+
+Matrix v2m(Vec3f v) {
+    Matrix m(4, 1);
+    m[0][0] = v.x;
+    m[1][0] = v.y;
+    m[2][0] = v.z;
+    m[3][0] = 1.f;
+    return m;
+}
+
+Matrix viewport(int x, int y, int w, int h) {
+    Matrix m = Matrix::identity(4);
+    m[0][3] = x+w/2.f;
+    m[1][3] = y+h/2.f;
+    m[2][3] = depth/2.f;
+
+    m[0][0] = w/2.f;
+    m[1][1] = h/2.f;
+    m[2][2] = depth/2.f;
+    return m;
+}
 
 void triangle(Vec3i t0, Vec3i t1, Vec3i t2, Vec2i uv0, Vec2i uv1, Vec2i uv2, TGAImage &image, float intensity, int *zbuffer) {
     if (t0.y==t1.y && t0.y==t2.y) return; // i dont care about degenerate triangles
@@ -57,6 +83,10 @@ int main(int argc, char** argv) {
     }
 
     { // draw the model
+        Matrix Projection = Matrix::identity(4);
+        Matrix ViewPort   = viewport(width/8, height/8, width*3/4, height*3/4);
+        Projection[3][2] = -1.f/camera.z;
+
         TGAImage image(width, height, TGAImage::RGB);
         for (int i=0; i<model->nfaces(); i++) {
             std::vector<int> face = model->face(i);
@@ -64,7 +94,7 @@ int main(int argc, char** argv) {
             Vec3f world_coords[3];
             for (int j=0; j<3; j++) {
                 Vec3f v = model->vert(face[j]);
-                screen_coords[j] = Vec3i((v.x+1.)*width/2., (v.y+1.)*height/2., (v.z+1.)*depth/2.);
+                screen_coords[j] =  m2v(ViewPort*Projection*v2m(v));
                 world_coords[j]  = v;
             }
             Vec3f n = (world_coords[2]-world_coords[0])^(world_coords[1]-world_coords[0]);