diff --git a/Makefile b/Makefile index 77bef81..d7b2900 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ TARGET = ray OBJECTS = main.o sampling.o camera.o environment.o disk.o plane.o sphere.o \ phong_brdf.o hsa_brdf.o directional_light.o point_light.o \ spot_light.o sphere_area_light.o disk_area_light.o scene.o tracer.o \ - path_tracer.o whitted_tracer.o + path_tracer.o whitted_tracer.o rgbe.o kd_tree.o DEPENDS = $(OBJECTS:.o=.d) CXXFLAGS = -ansi -pedantic -Wall -DGLM_FORCE_RADIANS -fopenmp LDLIBS = -lfreeimage -ljson_spirit diff --git a/kd_tree.cpp b/kd_tree.cpp new file mode 100644 index 0000000..23d544d --- /dev/null +++ b/kd_tree.cpp @@ -0,0 +1,313 @@ +#include "kd_tree.hpp" +#include + + +treeNode::treeNode(Photon p, superKey plane) +{ + photon = p; + haveChilds = false; + leftChild = NULL; + rightChild = NULL; + cutPlane = plane; +} + +treeNode::~treeNode() +{ + haveChilds = false; + if(leftChild) delete leftChild; + if(rightChild) delete rightChild; + leftChild = NULL; + rightChild = NULL; +} + +inline bool treeNode::operator==(treeNode b) +{ + return photon == b.photon; +} + +inline bool treeNode::operator==(Photon p) +{ + return photon == p; +} + +inline bool treeNode::operator>(treeNode b) +{ + return photon.greater(b.photon, cutPlane); +} + +inline bool treeNode::operator>(Photon p) +{ + return photon.greater(p, cutPlane); +} + +inline bool treeNode::operator<(treeNode b) +{ + return photon.lower(b.photon, cutPlane); +} + +inline bool treeNode::operator<(Photon p) +{ + return photon.lower(p, cutPlane); +} + +inline bool treeNode::operator>=(Vec3 p) +{ + return (cutPlane == XYZ && photon.position.x >= p.x) || + (cutPlane == YZX && photon.position.y >= p.y) || + (cutPlane == ZXY && photon.position.z >= p.z); +} + +inline bool treeNode::operator<=(Vec3 p) +{ + return (cutPlane == XYZ && photon.position.x <= p.x) || + (cutPlane == YZX && photon.position.y <= p.y) || + (cutPlane == ZXY && photon.position.z <= p.z); +} + +template +void copyArray(T a[], int size, T b[]) +{ + for(int i = 0; i < size; i++) + { + b[i] = a[i]; + } +} + +template +void restoreArray(T a, int size) +{ + for(int i = 0; i < size; i++) + a[i] = i; +} + +template +void Merge(T a[], int begin, int middle, int end, T b[], superKey key, std::vector originalData) +{ + int i = begin, j = middle; + + for(int k = begin; k < end; k++) + { + if(i < middle && (j >= end || originalData[a[i]].lessEqual(originalData[a[j]], key))) + { + b[k] = a[i]; + i++; + } + else + { + b[k] = a[j]; + j++; + } + } +} + +template +void SplitMerge(T b[], int begin, int end, T a[], superKey key, std::vector originalData) +{ + if(end - begin < 2) + return; + + int middle = (begin + end) / 2; + + SplitMerge(a, begin, middle, b, key, originalData); + SplitMerge(a, middle, end, b, key, originalData); + + Merge(b, begin, middle, end, a, key, originalData); +} + +template +void MegeSort(T a[], T b[], int size, superKey key, std::vector originalData) +{ + SplitMerge(b, 0, size, a, key ,originalData); +} + +kdTree::kdTree(){root = NULL;} +kdTree::~kdTree(){} + +void kdTree::addPhoton(Photon p) +{ + Photons.push_back(p); +} + + +void kdTree::createNodeKdTree(treeNode** node, std::vector originalData , int* xyz, int* yzx, int* zxy, superKey key, int begin, int end, int* xyz_2, int* yzx_2, int* zxy_2) +{ + if(end - begin < 2) + { + switch(key) + { + case XYZ: + *node = new treeNode(originalData[xyz[begin]], XYZ); + //std::cout << "Photon posicion " << xyz[begin] << " " << originalData[xyz[begin]] << std::endl; + break; + case YZX: + *node = new treeNode(originalData[yzx[begin]], YZX); + //std::cout << "Photon posicion " << yzx[begin] << " " << originalData[yzx[begin]] << std::endl; + break; + case ZXY: + *node = new treeNode(originalData[zxy[begin]], ZXY); + //std::cout << "Photon posicion " << zxy[begin] << " " << originalData[zxy[begin]] << std::endl; + break; + } + + return; + } + + int mid = (begin + end) / 2; + + switch(key) + { + case XYZ: + *node = new treeNode(originalData[xyz[mid]], XYZ); + (*node)->setChildsFlag(true); + //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); + key = YZX; + break; + case YZX: + *node = new treeNode(originalData[yzx[mid]], YZX); + (*node)->setChildsFlag(true); + //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); + key = ZXY; + break; + case ZXY: + *node = new treeNode(originalData[zxy[mid]], ZXY); + (*node)->setChildsFlag(true); + //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); + key = XYZ; + break; + } + + //std::cout<<"Rama izquierda" << std::endl; + createNodeKdTree((*node)->getLeftChild(), originalData, xyz_2, yzx_2, zxy_2, key, begin, mid, xyz, yzx, zxy); + //std::cout<<"Rama derecha" << std::endl; + createNodeKdTree((*node)->getRightChild(), originalData, xyz_2, yzx_2, zxy_2, key, mid + 1, end, xyz, yzx, zxy); +} + +void kdTree::reorderArrays(std::vector originalData, int* A1, int* A2, int begin, int mid, int end, int orderIndex, superKey key, int* B1, int* B2) +{ + int lowerindex1 = begin, higherindex1 = mid + 1, lowerindex2 = begin, higherindex2 = mid + 1; + + for(int i = begin; i < end; i++) + { + if(A1[i] != orderIndex) + { + if(originalData[A1[i]].lower(originalData[orderIndex], key)) + { + B1[lowerindex1] = A1[i]; + lowerindex1++; + } + else + { + B1[higherindex1] = A1[i]; + higherindex1++; + } + } + + if(A2[i] != orderIndex) + { + if(originalData[A2[i]].lower(originalData[orderIndex], key)) + { + B2[lowerindex2] = A2[i]; + lowerindex2++; + } + else + { + B2[higherindex2] = A2[i]; + higherindex2++; + } + } + } +} + +bool kdTree::buildKdTree() +{ + int size = Photons.size(); + //Arreglos con las superclaves de los photones + int *xyz = new int[size]; + int *yzx = new int[size]; + int *zxy = new int[size]; + int *xyz_aux = new int[size]; + int *yzx_aux = new int[size]; + int *zxy_aux = new int[size]; + + for(int i = 0; i < size; i++) + { + xyz[i] = i; + yzx[i] = i; + zxy[i] = i; + xyz_aux[i] = i; + } + + MegeSort(xyz, xyz_aux, size, XYZ, Photons); + restoreArray(xyz_aux, size); + + MegeSort(yzx, xyz_aux, size, YZX, Photons); + restoreArray(xyz_aux, size); + + MegeSort(zxy, xyz_aux, size, ZXY, Photons); + + for(int i = 0; i < size; i++) + { + xyz_aux[i] = xyz[i]; + yzx_aux[i] = yzx[i]; + zxy_aux[i] = zxy[i]; + } + + createNodeKdTree(&root, Photons , xyz, yzx, zxy, XYZ, 0, size, xyz_aux, yzx_aux, zxy_aux); + + //printTree(); + + return true; +} + +void kdTree::printTree(){ + std::queue q; + q.push(root); + + while(!q.empty()) + { + treeNode* n = q.front(); + q.pop(); + if(n) + { + std::cout << *n << " "; + + if(q.empty() || !q.front() || *(q.front()->getCutPlane()) != *(n->getCutPlane())) + std::cout << std::endl; + + q.push(*(n->getLeftChild())); + q.push(*(n->getRightChild())); + } + } +} + +void kdTree::printNode(treeNode* node) +{ +} + +std::vector kdTree::findInRange (Vec3 min, Vec3 max) +{ + std::vector photons; + + findInRange(min, max, photons, root); + + return photons; +} + +void kdTree::findInRange (Vec3 min, Vec3 max, std::vector &photons, treeNode *node) +{ + if(node == NULL) return; + + Vec3 position = node->getPhoton()->position; + if(position >= min && position <= max) + photons.push_back(*node->getPhoton()); + + if(*node >= min) + findInRange(min, max, photons, *node->getLeftChild()); + + if(*node <= max) + findInRange(min, max, photons, *node->getRightChild()); + +} diff --git a/kd_tree.hpp b/kd_tree.hpp new file mode 100644 index 0000000..89ef2d5 --- /dev/null +++ b/kd_tree.hpp @@ -0,0 +1,185 @@ +#include +#include +#include + +#include "rgbe.hpp" + +enum superKey + { + XYZ, + YZX, + ZXY + }; + +struct Vec3 +{ + float x; + float y; + float z; + + Vec3() { + x = y = z = 0.0f; + } + + inline bool equalFloat(const float x, const float y) + { + return (x - std::numeric_limits::epsilon() <= y) && (x + std::numeric_limits::epsilon() >= y); + } + + inline bool operator<=(const Vec3 b) + { + return (x < b.x || (equalFloat(x, b.x) && (y < b.y || (equalFloat(y, b.y) && z <= b.z)))); + } + + inline bool operator>=(const Vec3 b) + { + return (x > b.x || (equalFloat(x, b.x) && (y > b.y || (equalFloat(y, b.y) && z >= b.z)))); + } + + inline friend std::ostream& operator<<(std::ostream& out, const Vec3& v) + { + return out << "X:" << v.x << " Y:" << v.y << " Z:" << v.z; + } +}; + +struct Photon +{ + Vec3 position; + unsigned char radiance[4]; + char phi, theta; + short unused_flag; // for 20 bytes struct. + + Photon(Vec3 _p = Vec3(), float red = 0.0f, float green = 0.0f, float blue = 0.0f, char _phi = 0.0f, char _theta = 0.0f): + position(_p), + phi(_phi), + theta(_theta) + { + float2rgbe(radiance, red, green, blue); + } + + inline bool equalFloat(const float x, const float y) + { + return (x - std::numeric_limits::epsilon() <= y) && (x + std::numeric_limits::epsilon() >= y); + } + + inline bool operator==(const Photon p) + { + return equalFloat(position.x, p.position.x) && equalFloat(position.y, p.position.y) && equalFloat(position.z, p.position.z); + } + + inline bool lessEqual(const Photon& b, superKey key) + { + return (key == XYZ && (position.x < b.position.x || (equalFloat(position.x, b.position.x) && (position.y < b.position.y || (equalFloat(position.y, b.position.y) && position.z <= b.position.z))))) + || (key == YZX && (position.y < b.position.y || (equalFloat(position.y, b.position.y) && (position.z < b.position.z || (equalFloat(position.z, b.position.z) && position.x <= b.position.x))))) + || (key == ZXY && (position.z < b.position.z || (equalFloat(position.z, b.position.z) && (position.x < b.position.x || (equalFloat(position.x, b.position.x) && position.y <= b.position.y))))); + } + + inline bool lower(const Photon &b, superKey key) + { + return (key == XYZ && (position.x < b.position.x || (equalFloat(position.x, b.position.x) && (position.y < b.position.y || (equalFloat(position.y, b.position.y) && position.z < b.position.z))))) + || (key == YZX && (position.y < b.position.y || (equalFloat(position.y, b.position.y) && (position.z < b.position.z || (equalFloat(position.z, b.position.z) && position.x < b.position.x))))) + || (key == ZXY && (position.z < b.position.z || (equalFloat(position.z, b.position.z) && (position.x < b.position.x || (equalFloat(position.x, b.position.x) && position.y < b.position.y))))); + } + + inline bool greater(const Photon &b, superKey key) + { + return (key == XYZ && (position.x > b.position.x || (equalFloat(position.x, b.position.x) && (position.y > b.position.y || (equalFloat(position.y, b.position.y) && position.z > b.position.z))))) + || (key == YZX && (position.y > b.position.y || (equalFloat(position.y, b.position.y) && (position.z > b.position.z || (equalFloat(position.z, b.position.z) && position.x > b.position.x))))) + || (key == ZXY && (position.z > b.position.z || (equalFloat(position.z, b.position.z) && (position.x > b.position.x || (equalFloat(position.x, b.position.x) && position.y > b.position.y))))); + } + + inline friend std::ostream& operator<<(std::ostream& out, const Photon& p) + { + return out << "Position: " << p.position; + } +}; + +class treeNode +{ +private: + Photon photon; + bool haveChilds; + superKey cutPlane; + treeNode* leftChild; + treeNode* rightChild; +public: + treeNode(Photon p, superKey plane); + ~treeNode(); + inline bool operator==(treeNode b); + inline bool operator>(treeNode b); + inline bool operator<(treeNode b); + inline bool operator==(Photon p); + inline bool operator>(Photon p); + inline bool operator<(Photon p); + inline bool operator>=(Vec3 p); + inline bool operator<=(Vec3 p); + inline friend std::ostream& operator<<(std::ostream& out, treeNode &t) + { + out << "Cut Plane "; + + switch(*(t.getCutPlane())) + { + case XYZ: + out << "X "; + break; + case YZX: + out << "Y "; + break; + case ZXY: + out << "Z "; + break; + } + return out << "Photon " << *(t.getPhoton()) << " HaveChilds: " << *(t.getChildsFlag()); + } + + Photon* getPhoton(){return &photon;} + bool* getChildsFlag(){return &haveChilds;} + void setChildsFlag(bool childs){haveChilds = childs;} + superKey* getCutPlane(){return &cutPlane;} + treeNode** getLeftChild(){return &leftChild;} + treeNode** getRightChild(){return &rightChild;} +}; + +//Clase que genera el kdTree +class kdTree +{ +public: + kdTree(); + ~kdTree(); + + void addPhoton(Photon p); + bool buildKdTree(); + void printTree(); + std::vector findInRange (Vec3 min, Vec3 max); + +private: + treeNode* root; + std::vector Photons; + + void createNodeKdTree(treeNode** node, + std::vector originalData , + int* xyz, + int* yzx, + int* zxy, + superKey key, + int begin, + int end, + int* xyz_2, + int* yzx_2, + int* zxy_2); + + void reorderArrays(std::vector originalData, + int* A1, + int* A2, + int begin, + int mid, + int end, + int orderIndex, + superKey Key, + int* B1, + int* B2); + + void printNode(treeNode* node); + + void findInRange (Vec3 min, Vec3 max, std::vector &photons, treeNode *node); +}; diff --git a/rgbe.cpp b/rgbe.cpp new file mode 100644 index 0000000..c859ac5 --- /dev/null +++ b/rgbe.cpp @@ -0,0 +1,40 @@ +#include + +#include "rgbe.hpp" + +/* standard conversion from float pixels to rgbe pixels */ +/* note: you can remove the "inline"s if your compiler complains about it */ +void float2rgbe(unsigned char rgbe[4], float red, float green, float blue) { + float v; + int e; + + v = red; + if (green > v) v = green; + if (blue > v) v = blue; + if (v < 1e-32) { + rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0; + } + else { + v = frexp(v,&e) * 256.0/v; + rgbe[0] = (unsigned char) (red * v); + rgbe[1] = (unsigned char) (green * v); + rgbe[2] = (unsigned char) (blue * v); + rgbe[3] = (unsigned char) (e + 128); + } +} + +/* standard conversion from rgbe to float pixels */ +/* note: Ward uses ldexp(col+0.5,exp-(128+8)). However we wanted pixels */ +/* in the range [0,1] to map back into the range [0,1]. */ +void rgbe2float(float & red, float & green, float & blue, unsigned char rgbe[4]) { + float f; + + if (rgbe[3]) { /*nonzero pixel*/ + f = ldexp(1.0,rgbe[3]-(int)(128+8)); + red = rgbe[0] * f; + green = rgbe[1] * f; + blue = rgbe[2] * f; + } + else + red = green = blue = 0.0; +} diff --git a/rgbe.hpp b/rgbe.hpp new file mode 100644 index 0000000..5262b20 --- /dev/null +++ b/rgbe.hpp @@ -0,0 +1,27 @@ +/* THIS CODE CARRIES NO GUARANTEE OF USABILITY OR FITNESS FOR ANY PURPOSE. + * WHILE THE AUTHORS HAVE TRIED TO ENSURE THE PROGRAM WORKS CORRECTLY, + * IT IS STRICTLY USE AT YOUR OWN RISK. */ + +/* This file contains code to read and write four byte rgbe file format + developed by Greg Ward. It handles the conversions between rgbe and + pixels consisting of floats. The data is assumed to be an array of floats. + By default there are three floats per pixel in the order red, green, blue. + (RGBE_DATA_??? values control this.) Only the mimimal header reading and + writing is implemented. Each routine does error checking and will return + a status value as defined below. This code is intended as a skeleton so + feel free to modify it to suit your needs. + + (Place notice here if you modified the code.) + posted to http://www.graphics.cornell.edu/~bjw/ + written by Bruce Walter (bjw@graphics.cornell.edu) 5/26/95 + based on code written by Greg Ward +*/ + +#pragma once +#ifndef RGBE_HPP +#define RGBE_HPP + +extern void float2rgbe(unsigned char rgbe[4], float red, float green, float blue); +extern void rgbe2float(float & red, float & green, float & blue, unsigned char rgbe[4]); + +#endif