diff --git a/kd_tree.cpp b/kd_tree.cpp
index affad6a..e733d14 100644
--- a/kd_tree.cpp
+++ b/kd_tree.cpp
@@ -1,17 +1,21 @@
-#ifndef NDEBUG
 #include <iostream>
+#ifndef NDEBUG
 #include <fstream>
 #endif
 #include <queue>
 
+#include <omp.h>
+
 #include "kd_tree.hpp"
 
 #ifndef NDEBUG
 using std::ofstream;
 using std::ios;
-using std::endl;
 #endif
 
+using std::cout;
+using std::endl;
+
 treeNode::treeNode(Photon p, superKey plane)
 {
   photon = p;
@@ -86,12 +90,13 @@ void copyArray(T a[], int size, T b[])
 template <class T>
 void restoreArray(T a, int size)
 {
+#pragma omp parallel for schedule(dynamic, 1)
   for(int i = 0; i < size; i++)
     a[i] = i;
 }
 
 template <class T>
-void Merge(T a[], int begin, int middle, int end, T b[], superKey key, std::vector<Photon> originalData)
+void Merge(T a[], int begin, int middle, int end, T b[], superKey key, std::vector<Photon> & originalData)
 {
   int i = begin, j = middle;
 
@@ -111,7 +116,7 @@ void Merge(T a[], int begin, int middle, int end, T b[], superKey key, std::vect
 }
 
 template <class T>
-void SplitMerge(T b[], int begin, int end, T a[], superKey key, std::vector<Photon> originalData)
+void SplitMerge(T b[], int begin, int end, T a[], superKey key, std::vector<Photon> & originalData)
 {
   if(end - begin < 2)
     return;
@@ -125,9 +130,9 @@ void SplitMerge(T b[], int begin, int end, T a[], superKey key, std::vector<Phot
 }
 
 template <class T>
-void MegeSort(T a[], T b[], int size, superKey key, std::vector<Photon> originalData)
+void MegeSort(T a[], T b[], int size, superKey key, std::vector<Photon> & originalData)
 {
-  SplitMerge(b, 0, size, a, key ,originalData);
+  SplitMerge(b, 0, size, a, key, originalData);
 }
 
 kdTree::kdTree(){root = NULL;}
@@ -240,25 +245,87 @@ bool kdTree::buildKdTree()
   int *yzx = new int[size];
   int *zxy = new int[size];
   int *xyz_aux = new int[size];
+  int *xyz_aux2 = new int[size];
+  int *xyz_aux3 = new int[size];
   int *yzx_aux = new int[size];
   int *zxy_aux = new int[size];
 
+  cout << "Calculating medians." << endl;
+
+#pragma omp parallel for schedule(dynamic, 1)
   for(int i = 0; i < size; i++)
     {
       xyz[i] = i;
       yzx[i] = i;
       zxy[i] = i;
       xyz_aux[i] = i;
+      xyz_aux2[i] = i;
+      xyz_aux3[i] = i;
     }
 
-  MegeSort(xyz, xyz_aux, size, XYZ, Photons);
-  restoreArray(xyz_aux, size);
+  if (omp_get_max_threads() == 2) {
 
-  MegeSort(yzx, xyz_aux, size, YZX, Photons);
-  restoreArray(xyz_aux, size);
+#pragma omp parallel
+    {
+      if (omp_get_thread_num() == 0) {
+#pragma omp critical
+	{
+	  cout << "Sorting \x1b[1;33mXYZ\x1b[m." << endl;
+	}
+	MegeSort(xyz, xyz_aux, size, XYZ, Photons);
+	// restoreArray(xyz_aux, size);
+      } else if(omp_get_thread_num() == 1) {
+#pragma omp critical
+	{
+	  cout << "Sorting \x1b[1;33mYZX\x1b[m." << endl;
+	}
+	MegeSort(yzx, xyz_aux2, size, YZX, Photons);
+	// restoreArray(xyz_aux, size);
+      }
+    }
 
-  MegeSort(zxy, xyz_aux, size, ZXY, Photons);
+    cout << "Sorting \x1b[1;33mZXY\x1b[m." << endl;
+    MegeSort(zxy, xyz_aux3, size, ZXY, Photons);
+    
+  } else if (omp_get_max_threads() >= 3) {
+#pragma omp parallel
+    {
+      if (omp_get_thread_num() == 0) {
+#pragma omp critical
+	{
+	  cout << "Sorting \x1b[1;33mXYZ\x1b[m." << endl;
+	}
+	MegeSort(xyz, xyz_aux, size, XYZ, Photons);
+	// restoreArray(xyz_aux, size);
+      } else if(omp_get_thread_num() == 1) {
+#pragma omp critical
+	{
+	  cout << "Sorting \x1b[1;33mYZX\x1b[m." << endl;
+	}
+	MegeSort(yzx, xyz_aux2, size, YZX, Photons);
+	// restoreArray(xyz_aux, size);
+      } else if (omp_get_thread_num() == 2) {
+#pragma omp critical
+	{
+	  cout << "Sorting \x1b[1;33mZXY\x1b[m." << endl;
+	}
+	MegeSort(zxy, xyz_aux3, size, ZXY, Photons);
+      }
+    }
+  } else {
+    cout << "Sorting \x1b[1;33mXYZ\x1b[m." << endl;
+    MegeSort(xyz, xyz_aux, size, XYZ, Photons);
+    // restoreArray(xyz_aux, size);
+    
+    cout << "Sorting \x1b[1;33mYZX\x1b[m." << endl;
+    MegeSort(yzx, xyz_aux2, size, YZX, Photons);
+    // restoreArray(xyz_aux, size);
+    
+    cout << "Sorting \x1b[1;33mZXY\x1b[m." << endl;
+    MegeSort(zxy, xyz_aux3, size, ZXY, Photons);
+  }
 
+#pragma omp parallel for schedule(dynamic, 1)
   for(int i = 0; i < size; i++)
     {
       xyz_aux[i] = xyz[i];
@@ -266,11 +333,13 @@ bool kdTree::buildKdTree()
       zxy_aux[i] = zxy[i];
     }
 
+  cout << "Adding photons to the tree." << endl;
   createNodeKdTree(&root, Photons , xyz, yzx, zxy, XYZ, 0, size, xyz_aux, yzx_aux, zxy_aux);
 
   //printTree();
 
 #ifndef NDEBUG
+  cout << "Writing photons to \x1b[1;33mphotons.txt\x1b[m" << endl;
   ofstream ofs("photons.txt", ios::out);
   float r, g, b;
   for (std::vector<Photon>::iterator it = Photons.begin(); it != Photons.end(); it++) {
@@ -279,6 +348,15 @@ bool kdTree::buildKdTree()
   }
   ofs.close();
 #endif
+
+  // delete[] xyz;
+  // delete[] yzx;
+  // delete[] zxy;
+  // delete[] xyz_aux;
+  // delete[] xyz_aux2;
+  // delete[] xyz_aux3;
+  // delete[] yzx_aux;
+  // delete[] zxy_aux;
   
   return true;
 }
@@ -332,3 +410,7 @@ void kdTree::findInRange (Vec3 min, Vec3 max, std::vector<Photon> &photons, tree
     findInRange(min, max, photons, *node->getRightChild());
 
 }
+
+size_t kdTree::getNumPhotons() {
+  return Photons.size();
+}
diff --git a/kd_tree.hpp b/kd_tree.hpp
index 3465cd7..d75fc14 100644
--- a/kd_tree.hpp
+++ b/kd_tree.hpp
@@ -154,6 +154,8 @@ public:
   void printTree();
   std::vector<Photon> findInRange (Vec3 min, Vec3 max) const;
 
+  size_t getNumPhotons();
+  
 private:
   treeNode* root;
   std::vector<Photon> Photons;
diff --git a/main.cpp b/main.cpp
index 10c5e31..8a0ec18 100644
--- a/main.cpp
+++ b/main.cpp
@@ -27,6 +27,8 @@ using namespace glm;
 ////////////////////////////////////////////
 #define ANSI_BOLD_YELLOW "\x1b[1;33m"
 #define ANSI_RESET_STYLE "\x1b[m"
+#define MAX_W 1920
+#define MAX_H 1080
 
 ////////////////////////////////////////////
 // Function prototypes.
@@ -51,7 +53,7 @@ static float g_fov = 45.0f;
 static int g_w = 640;
 static int g_h = 480;
 static float g_a_ratio = 640.0f / 480.0f;
-static vec3 ** image;
+static vec3 image[MAX_H][MAX_W];
 static tracer_t g_tracer = NONE;
 static unsigned int g_max_depth = 5;
 static float g_gamma = 2.2f;
@@ -82,11 +84,6 @@ int main(int argc, char ** argv) {
   // Initialize everything.
   FreeImage_Initialise();
 
-  image = new vec3*[g_h];
-  for (int i = 0; i < g_h; i++) {
-    image[i] = new vec3[g_w];
-  }
-
   try {
     scn = new Scene(g_input_file);
   } catch (SceneError & e) {
@@ -198,10 +195,6 @@ int main(int argc, char ** argv) {
   delete scn;
   delete tracer;
 
-  for (int i = 0; i < g_h; i++)
-    delete[] image[i];
-  delete[] image;
-
   FreeImage_DeInitialise();
   
   return EXIT_SUCCESS;
@@ -227,6 +220,8 @@ void print_usage(char ** const argv) {
   cerr << "    \tDefaults to 25 samples." << endl;
   cerr << "  -w\tImage size in pixels as \"WIDTHxHEIGHT\"." << endl;
   cerr << "    \tDefaults to 640x480 pixels." << endl;
+  cerr << "    \tMinimum resolution is 1x1 pixels." << endl;
+  cerr << "    \tMaxmimum resolution is " << MAX_W << "x" << MAX_H << " pixels." << endl;
   cerr << "  -r\tMaxmimum recursion depth." << endl;
   cerr << "    \tDefaults to 5." << endl;
   cerr << "  -g\tGamma correction value (>= 0)." << endl;
@@ -295,7 +290,7 @@ void parse_args(int argc, char ** const argv) {
 	optarg[x_pos] = '\0';
 	g_w = atoi(optarg);
 	g_h = atoi(&optarg[x_pos + 1]);
-	if (g_w <= 0 || g_h <= 0) {
+	if (g_w <= 0 || g_h <= 0 || g_w >= MAX_W || g_h >= MAX_H) {
 	  cerr << "Invalid screen resolution: " << optarg << endl;
 	  print_usage(argv);
 	  exit(EXIT_FAILURE);
diff --git a/photon_tracer.cpp b/photon_tracer.cpp
index 563338c..bd840b6 100644
--- a/photon_tracer.cpp
+++ b/photon_tracer.cpp
@@ -115,8 +115,8 @@ vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
 	p_contrib += ((_f->m_mat->m_diffuse / pi<float>()) * vec3(red, green, blue)) / (pi<float>() * (radius * radius));
       }
       
-      color += (1.0f - _f->m_mat->m_rho) * (((dir_diff_color + p_contrib) * (_f->m_mat->m_diffuse / pi<float>())) +
-					    (_f->m_mat->m_specular * dir_spec_color));
+      color += (1.0f - _f->m_mat->m_rho) * (((dir_diff_color) * (_f->m_mat->m_diffuse / pi<float>())) +
+					    (_f->m_mat->m_specular * dir_spec_color) + p_contrib);
 
       // Determine the specular reflection color.
       if (_f->m_mat->m_rho > 0.0f && rec_level < m_max_depth) {
@@ -211,6 +211,7 @@ void PhotonTracer::build_photon_map(Scene * s, const size_t n_photons_per_ligth,
   }
   cout << endl;
 
+  cout << "Generated " << ANSI_BOLD_YELLOW << m_photon_map.getNumPhotons() << ANSI_RESET_STYLE << " total photons." << endl;
   cout << "Building photon map Kd-tree." << endl;
   m_photon_map.buildKdTree();
 }