Added "fixes" of the kd tree.

Makefile

@@ -6,7 +6,7 @@ OBJECTS = main.o sampling.o camera.o environment.o disk.o plane.o sphere.o \
           spot_light.o sphere_area_light.o disk_area_light.o scene.o tracer.o \
           path_tracer.o whitted_tracer.o rgbe.o kd_tree.o photon_tracer.o
 DEPENDS = $(OBJECTS:.o=.d)
-CXXFLAGS = -ansi -pedantic -Wall -DGLM_FORCE_RADIANS -fopenmp -std=c++11
+CXXFLAGS = -std=c++11 -pedantic -Wall -DGLM_FORCE_RADIANS -fopenmp -D_USE_CPP11_RANDOM #-DENABLE_KD_TREE
 LDLIBS = -lfreeimage -ljson_spirit
 
 .PHONY: all

disk_area_light.hpp

@@ -7,7 +7,9 @@
 
 class DiskAreaLight: public AreaLight {
 public:
-  DiskAreaLight(Disk * _s, float _c = 1.0, float _l = 0.0, float _q = 0.0): AreaLight(static_cast<Figure *>(_s), _c, _l, _q) { }
+  DiskAreaLight(Disk * _s, float _c = 1.0, float _l = 0.0, float _q = 0.0):
+    AreaLight(static_cast<Figure *>(_s), _c, _l, _q)
+  { }
 
   virtual vec3 sample_at_surface();
 };

kd_tree.cpp

@@ -139,25 +139,66 @@ void kdTree::addPhoton(Photon p)
   Photons.push_back(p);
 }
 
-void kdTree::createNodeKdTree(treeNode** node, std::vector<Photon> & originalData , int* xyz, int* yzx, int* zxy, superKey key, int begin, int end, int* xyz_2, int* yzx_2, int* zxy_2)
+void kdTree::createNodeKdTree(treeNode** node, std::vector<Photon> originalData , int xyz[], int yzx[], int zxy[], superKey key, int begin, int end, int fullsize)
 {
-  if(end - begin < 2)
+
+  int size = end - begin;
+
+  if(size <= 2)
+  {
+    if(size < 2)
     {
       switch(key)
       {
         case XYZ:
           *node = new treeNode(originalData[xyz[begin]], XYZ);
-	  //std::cout << "Photon posicion " << xyz[begin] << " " << originalData[xyz[begin]] << std::endl;
+          if(size == 1)
+          {
+            (*node)->setChildsFlag(true);
+            *((*node)->getRightChild()) = new treeNode(originalData[xyz[end]], YZX);
+          }
         break;
         case YZX:
           *node = new treeNode(originalData[yzx[begin]], YZX);
-	  //std::cout << "Photon posicion " << yzx[begin] << " " << originalData[yzx[begin]] << std::endl;
+          if(size == 1)
+          {
+            (*node)->setChildsFlag(true);
+            *((*node)->getRightChild()) = new treeNode(originalData[yzx[end]], ZXY);
+          }
         break;
         case ZXY:
           *node = new treeNode(originalData[zxy[begin]], ZXY);
-	  //std::cout << "Photon posicion " << zxy[begin] << " " << originalData[zxy[begin]] << std::endl;
+          if(size == 1)
+          {
+            (*node)->setChildsFlag(true);
+            *((*node)->getRightChild()) = new treeNode(originalData[zxy[end]], XYZ);
+          }
         break;
       }
+    }
+    else
+    {
+      int mid = (begin + end) / 2;
+
+      switch(key)
+      {
+        case XYZ:
+          *node = new treeNode(originalData[xyz[mid]], XYZ);
+          *((*node)->getLeftChild()) = new treeNode(originalData[xyz[begin]], YZX);
+          *((*node)->getRightChild()) = new treeNode(originalData[xyz[end]], YZX);
+        break;
+        case YZX:
+          *node = new treeNode(originalData[yzx[mid]], YZX);
+          *((*node)->getLeftChild()) = new treeNode(originalData[yzx[begin]], ZXY);
+          *((*node)->getRightChild()) = new treeNode(originalData[yzx[end]], ZXY);
+        break;
+        case ZXY:
+          *node = new treeNode(originalData[zxy[mid]], ZXY);
+          *((*node)->getLeftChild()) = new treeNode(originalData[zxy[begin]], XYZ);
+          *((*node)->getRightChild()) = new treeNode(originalData[zxy[end]], XYZ);
+        break;
+      }
+    }
 
     return;
   }
@@ -169,66 +210,135 @@ void kdTree::createNodeKdTree(treeNode** node, std::vector<Photon> & originalDat
     case XYZ:
       *node = new treeNode(originalData[xyz[mid]], XYZ);
       (*node)->setChildsFlag(true);
+      //std::cout << "CUT XYZ: " << xyz[mid] << std::endl;
       //std::cout << "Photon posicion " << xyz[mid] << " " << originalData[xyz[mid]] << std::endl;
-      reorderArrays(originalData, yzx, zxy, begin, mid, end, xyz[mid], XYZ, yzx_2, zxy_2);
+      reorderArrays(originalData, yzx, zxy, begin, mid, end, xyz[mid], XYZ);
+      xyz[mid] = -1;
       key = YZX;
       break;
     case YZX:
       *node = new treeNode(originalData[yzx[mid]], YZX);
       (*node)->setChildsFlag(true);
+      //std::cout << "CUT YZX: " << yzx[mid] << std::endl;
       //std::cout << "Photon posicion " << yzx[mid] << " " << originalData[yzx[mid]] << std::endl;
-      reorderArrays(originalData, xyz, zxy, begin, mid, end, yzx[mid], YZX, xyz_2, zxy_2);
+      reorderArrays(originalData, xyz, zxy, begin, mid, end, yzx[mid], YZX);
+      yzx[mid] = -1;
       key = ZXY;
       break;
     case ZXY:
       *node = new treeNode(originalData[zxy[mid]], ZXY);
       (*node)->setChildsFlag(true);
+      //std::cout << "CUT ZXY: " << zxy[mid] << std::endl;
       //std::cout << "Photon posicion " << zxy[mid] << " " << originalData[zxy[mid]] << std::endl;
-      reorderArrays(originalData, xyz, yzx, begin, mid, end, zxy[mid], ZXY, xyz_2, yzx_2);
+      reorderArrays(originalData, xyz, yzx, begin, mid, end, zxy[mid], ZXY);
+      zxy[mid] = -1;
       key = XYZ;
       break;
   }
 
+  /*std::cout << "XYZ: ";
+  for(int i = 0; i < fullsize; i++)
+    std::cout << xyz[i]  << " ";
+
+  std::cout << std::endl << "YZX: ";
+  for(int i = 0; i < fullsize; i++)
+    std::cout << yzx[i]  << " ";
+  
+  std::cout << std::endl << "ZXY: ";
+  for(int i = 0; i < fullsize; i++)
+    std::cout << zxy[i]  << " ";
+
+  std::cout << std::endl << std::endl;*/
+
   //std::cout<<"Rama izquierda" << std::endl;
-  createNodeKdTree((*node)->getLeftChild(), originalData, xyz_2, yzx_2, zxy_2, key, begin, mid, xyz, yzx, zxy);
+  createNodeKdTree((*node)->getLeftChild(), originalData, xyz, yzx, zxy, key, begin, mid - 1, fullsize);
   //std::cout<<"Rama derecha" << std::endl;
-  createNodeKdTree((*node)->getRightChild(), originalData, xyz_2, yzx_2, zxy_2, key, mid + 1, end, xyz, yzx, zxy);
+  createNodeKdTree((*node)->getRightChild(), originalData, xyz, yzx, zxy, key, mid + 1, end, fullsize);
 }
 
-void kdTree::reorderArrays(std::vector<Photon> & originalData, int* A1, int* A2, int begin, int mid, int end, int orderIndex, superKey key, int* B1, int* B2)
+void kdTree::reorderArrays(std::vector<Photon> originalData, int A1[], int A2[], int begin, int mid, int end, int orderIndex, superKey key)
 {
-  int lowerindex1 = begin, higherindex1 = mid + 1, lowerindex2 = begin, higherindex2 = mid + 1;
+  std::vector<int> lower1, lower2, higher1, higher2;
 
-  for(int i = begin; i < end; i++)
+  for(int i = begin; i <= end; i++)
   {   
     if(A1[i] != orderIndex)
     {
-	  if(originalData[A1[i]].lower(originalData[orderIndex], key))
+      if(originalData[A1[i]].lessEqual(originalData[orderIndex], key))
-	    {
+        lower1.push_back(A1[i]);
-	      B1[lowerindex1] = A1[i];
-	      lowerindex1++;
-	    }
       else
-	    {
+        higher1.push_back(A1[i]);
-	      B1[higherindex1] = A1[i];
-	      higherindex1++;
-	    }
     }
 
     if(A2[i] != orderIndex)
     {
-	  if(originalData[A2[i]].lower(originalData[orderIndex], key))
+      if(originalData[A2[i]].lessEqual(originalData[orderIndex], key))
-	    {
+        lower2.push_back(A2[i]);
-	      B2[lowerindex2] = A2[i];
-	      lowerindex2++;
-	    }
       else
+        higher2.push_back(A2[i]);
+    }
+  }
+
+  std::copy(lower1.begin(), lower1.end(), &A1[begin]);
+  std::copy(higher1.begin(), higher1.end(), &A1[mid + 1]);
+
+  std::copy(lower2.begin(), lower2.end(), &A2[begin]);
+  std::copy(higher2.begin(), higher2.end(), &A2[mid + 1]);
+  
+  A1[mid] = A2[mid] = -1;
+
+  lower1.clear();
+  lower2.clear();
+  higher1.clear();
+  higher2.clear();
+}
+
+void kdTree::removeDuplicates(std::vector<Photon> &originalData, int* &xyz, int* &yzx, int* &zxy, int &size)
 {
-	      B2[higherindex2] = A2[i];
+  if(size == 0) return;
-	      higherindex2++;
+
-	    }
+  std::vector<int> xyzNew;
+  std::vector<int> yzxNew;
+  std::vector<int> zxyNew;
+
+  xyzNew.push_back(xyz[0]);
+  yzxNew.push_back(yzx[0]);
+  zxyNew.push_back(zxy[0]);
+
+  for(int i = 1; i < size; i++)
+  {
+    if( !(originalData[xyz[i]] == originalData[xyz[i - 1]]) )
+      xyzNew.push_back(xyz[i]);
+    
+    if( !(originalData[yzx[i]] == originalData[yzx[i - 1]]) )
+      yzxNew.push_back(yzx[i]);
+    
+    if( !(originalData[zxy[i]] == originalData[zxy[i - 1]]) )
+      zxyNew.push_back(zxy[i]);
   }
+
+  if(xyzNew.size() != size)
+  {
+    std::cout << size - xyzNew.size() << " duplicates removed" << std::endl;
+
+    delete []xyz;
+    delete []yzx;
+    delete []zxy;
+
+    size = xyzNew.size();
+
+    xyz = new int[size];
+    yzx = new int[size];
+    zxy = new int[size];
+
+    std::copy(xyzNew.begin(), xyzNew.end(), &xyz[0]);
+    std::copy(yzxNew.begin(), yzxNew.end(), &yzx[0]);
+    std::copy(zxyNew.begin(), zxyNew.end(), &zxy[0]);
   }
+
+  xyzNew.clear();
+  yzxNew.clear();
+  zxyNew.clear();
 }
 
 bool kdTree::buildKdTree()
@@ -238,108 +348,36 @@ bool kdTree::buildKdTree()
   int *xyz = new int[size];
   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;
+  int *aux = new int[size];
 
-#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;
+    aux[i] = i;
-      xyz_aux2[i] = i;
-      xyz_aux3[i] = i;
   }
 
-  if (omp_get_max_threads() == 2) {
+  std::cout << "Initial sorting of data" << std::endl;
 
-#pragma omp parallel
+  MegeSort(xyz, aux, size, XYZ, Photons);
-    {
+  restoreArray(aux, size);
-      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);
-      }
-    }
 
-    cout << "Sorting \x1b[1;33mZXY\x1b[m." << endl;
+  MegeSort(yzx, aux, size, YZX, Photons);
-    MegeSort(zxy, xyz_aux3, size, ZXY, Photons);
+  restoreArray(aux, size);
 
-  } else if (omp_get_max_threads() >= 3) {
+  MegeSort(zxy, aux, size, ZXY, Photons);
-#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;
+  std::cout << "Removing Doubles" << std::endl;
-    MegeSort(yzx, xyz_aux2, size, YZX, Photons);
+  removeDuplicates(Photons, xyz, yzx, zxy, size);
-    // restoreArray(xyz_aux, size);
 
-    cout << "Sorting \x1b[1;33mZXY\x1b[m." << endl;
+  createNodeKdTree(&root, Photons , xyz, yzx, zxy, XYZ, 0, size - 1, size);
-    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];
-      yzx_aux[i] = yzx[i];
-      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();
 
   delete []xyz;
   delete []yzx;
   delete []zxy;
-  delete[] xyz_aux;
+  delete []aux;
-  delete[] xyz_aux2;
-  delete[] xyz_aux3;
-  delete[] yzx_aux;
-  delete[] zxy_aux;
  
   return true;
 }

kd_tree.hpp

@@ -9,7 +9,7 @@ enum superKey
   {
     XYZ,
     YZX,
-    ZXY
+    ZXY,
   };
 
 struct Vec3
@@ -172,28 +172,9 @@ private:
   treeNode* root;
   std::vector<Photon> Photons;
 
-  void createNodeKdTree(treeNode** node,
+  void createNodeKdTree(treeNode** node, std::vector<Photon> originalData , int xyz[], int yzx[], int zxy[], superKey key, int begin, int end, int fullsize);
-			std::vector<Photon> & originalData ,
+  void reorderArrays(std::vector<Photon> originalData, int A1[], int A2[], int begin, int mid, int end, int orderIndex, superKey key);
-			int* xyz,
+  void removeDuplicates(std::vector<Photon> &originalData, int* &xyz, int* &yzx, int* &zxy, int &size);
-			int* yzx,
-			int* zxy,
-			superKey key,
-			int begin,
-			int end,
-			int* xyz_2,
-			int* yzx_2,
-			int* zxy_2);
-
-  void reorderArrays(std::vector<Photon> & originalData,
-		     int* A1,
-		     int* A2,
-		     int begin,
-		     int mid,
-		     int end,
-		     int orderIndex,
-		     superKey Key,
-		     int* B1,
-		     int* B2);
 
   void printNode(treeNode* node);
 

photon_tracer.cpp

@@ -111,14 +111,24 @@ vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
 
       // Calculate photon map contribution
       radius = m_h_radius;
+#ifdef ENABLE_KD_TREE
+      Vec3 vmin(i_pos.x - m_h_radius, i_pos.y - m_h_radius, i_pos.z - m_h_radius);
+      Vec3 vmax(i_pos.x + m_h_radius, i_pos.y + m_h_radius, i_pos.z + m_h_radius);
+      photons = m_photon_map.findInRange(vmin, vmax);
+#else
       m_photon_map.find_by_distance(photons, i_pos, n, m_h_radius, 1000);
+#endif
       while(photons.size() == 0 && radius < 5.0) {
 	radius *= 2;
 	m_photon_map.find_by_distance(photons, i_pos, n, m_h_radius, 1000);
       }
 
       radius = m_h_radius;
+#ifdef ENABLE_KD_TREE
+      caustics = m_photon_map.findInRange(vmin, vmax);
+#else
       m_caustics_map.find_by_distance(caustics, i_pos, n, m_h_radius, 1000);
+#endif
       while(caustics.size() == 0 && radius < 5.0) {
 	radius *= 2;
 	m_caustics_map.find_by_distance(caustics, i_pos, n, m_h_radius, 1000);
@@ -237,11 +247,14 @@ void PhotonTracer::photon_tracing(Scene * s, const size_t n_photons_per_ligth, c
 
     assert(pl != NULL || al != NULL);
     
-#pragma omp parallel for schedule(dynamic, 1) private(l_sample, s_normal, h_sample, r1, r2) shared(current)
+#pragma omp parallel for schedule(dynamic, 1) private(l_sample, s_normal, h_sample, r1, r2, power, ls, dir, ph) shared(al, pl, current)
     for (size_t p = 0; p < n_photons_per_ligth; p++) {
       if (al != NULL) {
+#pragma omp critical
+	{
 	  l_sample = al->sample_at_surface();
 	  s_normal = al->normal_at_last_sample();
+	}
 	l_sample = l_sample + (BIAS * s_normal);
 	
 	if (!specular || (specular && spec_figures.size() == 0)) {
@@ -319,10 +332,12 @@ void PhotonTracer::build_photon_map(const char * photons_file, const bool causti
 
 void PhotonTracer::build_photon_map(const bool caustics) {
   cout << "Building photon map Kd-tree." << endl;
+#ifdef ENABLE_KD_TREE
   if (!caustics)
     m_photon_map.buildKdTree();
   else
   m_caustics_map.buildKdTree();
+#endif
 }
 
 void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_level) {

sampling.cpp

@@ -1,15 +1,23 @@
+#ifdef _USE_CPP11_RANDOM
+#warning "Using C++11 random number generators."
 #include <random>
 #include <chrono>
 #include <functional>
+#else
+#include <cstdlib>
+#include <ctime>
+#endif
 
 #include <glm/glm.hpp>
 #include <glm/gtc/constants.hpp>
 
 #include "sampling.hpp"
 
+#ifdef _USE_CPP11_RANDOM
 using std::uniform_real_distribution;
 using std::mt19937;
 using std::bind;
+#endif
 using glm::mat3;
 using glm::abs;
 using glm::normalize;
@@ -20,14 +28,27 @@ using glm::pi;
 const float PDF = (1.0f / (2.0f * pi<float>()));
 
 static bool seeded = false;
+
+#ifdef _USE_CPP11_RANDOM
 static uniform_real_distribution<float> dist(0, 1);
 static mt19937 engine;
 static auto generator = bind(dist, engine);
+#endif
 
 float random01() {
-  if (!seeded)
+  if (!seeded) {
+#ifdef _USE_CPP11_RANDOM
     engine.seed(std::chrono::system_clock::now().time_since_epoch().count());
+#else
+    srand(time(NULL));
+#endif
+    seeded = true;
+  }
+#ifdef _USE_CPP11_RANDOM
   return generator();
+#else
+  return static_cast<float>(rand()) / RAND_MAX;
+#endif
 }
 
 vec2 sample_pixel(int i, int j, float w, float h, float a_ratio, float fov) {