# HG changeset patch # User Radek Brich # Date 1208939913 -7200 # Node ID 9569e9f35374bbaafc63001f6f7c83418aab02ba # Parent dbe8438d5dcadd0f6784d8ad5c5bf4942c96349a move shapes to extra source file add serialize header and source file with common serialization functions dump/load feature for shapes and kd-tree fix few minor bugs diff -r dbe8438d5dca -r 9569e9f35374 TODO --- a/TODO Tue Apr 22 13:33:12 2008 +0200 +++ b/TODO Wed Apr 23 10:38:33 2008 +0200 @@ -7,10 +7,6 @@ * textures (3D procedural, pixmaps) * generalization: Camera "shader" (ray generator), surface shader and maybe light & background shaders * namespace - * kd-tree: - - optimize structures - - optimize construction: use box-shape intersection instead of bounding boxes of shapes - - save/load * Python binding for all classes * stochastic oversampling * absorbtion of refracted rays in dense materials (can be computed using shape distance and some 'absorbance' constant) diff -r dbe8438d5dca -r 9569e9f35374 include/common.h --- a/include/common.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/common.h Wed Apr 23 10:38:33 2008 +0200 @@ -31,6 +31,9 @@ #include #include #include +#include + +using namespace std; #ifdef PYRIT_DOUBLE # define Float double @@ -107,4 +110,19 @@ } } +inline void trim(string& s) +{ + const char *ws = " \n"; + string::size_type pos = s.find_last_not_of(ws); + if (pos != string::npos) + { + s.erase(pos + 1); + pos = s.find_first_not_of(ws); + if (pos != string::npos) + s.erase(0, pos); + } + else + s.erase(s.begin(), s.end()); +} + #endif diff -r dbe8438d5dca -r 9569e9f35374 include/container.h --- a/include/container.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/container.h Wed Apr 23 10:38:33 2008 +0200 @@ -29,7 +29,7 @@ #include -#include "scene.h" +#include "shapes.h" using namespace std; @@ -49,6 +49,10 @@ virtual Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray, Float &nearest_distance); virtual void optimize() {}; + + ShapeList & getShapes() { return shapes; }; + + virtual ostream & dump(ostream &st); }; #endif diff -r dbe8438d5dca -r 9569e9f35374 include/kdtree.h --- a/include/kdtree.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/kdtree.h Wed Apr 23 10:38:33 2008 +0200 @@ -53,19 +53,19 @@ ~KdNode(); void setLeaf() { flags |= 3; }; - bool isLeaf() { return (flags & 3) == 3; }; + const bool isLeaf() const { return (flags & 3) == 3; }; void setAxis(int aAxis) { flags &= ~3; flags |= aAxis; }; - short getAxis() { return flags & 3; }; + const int getAxis() const { return flags & 3; }; void setSplit(Float aSplit) { split = aSplit; }; - Float& getSplit() { return split; }; + const Float& getSplit() const { return split; }; void setChildren(KdNode *node) { children = node; assert((flags & 3) == 0); }; - KdNode* getLeftChild() { return (KdNode*)((off_t)children & ~3); }; - KdNode* getRightChild() { return (KdNode*)((off_t)children & ~3) + 1; }; + KdNode* getLeftChild() const { return (KdNode*)((off_t)children & ~3); }; + KdNode* getRightChild() const { return (KdNode*)((off_t)children & ~3) + 1; }; - ShapeList* getShapes() { return (ShapeList*)((off_t)shapes & ~3); }; + ShapeList* getShapes() const { return (ShapeList*)((off_t)shapes & ~3); }; void addShape(Shape* aShape) { getShapes()->push_back(aShape); }; }; @@ -79,6 +79,7 @@ int max_depth; void recursive_build(KdNode *node, BBox bbox, int maxdepth); + void recursive_load(istream &st, KdNode *node); public: KdTree() : Container(), root(NULL), built(false), max_depth(32) {}; ~KdTree() { if (root) delete root; }; @@ -87,9 +88,12 @@ Float &nearest_distance); void optimize() { build(); }; void build(); - void save(ostream &str, KdNode *node = NULL); - void load(istream &str, KdNode *node = NULL); + const bool isBuilt() const { return built; }; + KdNode *getRootNode() const { return root; }; void setMaxDepth(int md) { max_depth = md; }; + + ostream & dump(ostream &st); + istream & load(istream &st); }; #endif diff -r dbe8438d5dca -r 9569e9f35374 include/raytracer.h --- a/include/raytracer.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/raytracer.h Wed Apr 23 10:38:33 2008 +0200 @@ -34,8 +34,6 @@ #include "container.h" #include "scene.h" -using namespace std; - class Raytracer; struct RenderrowData { diff -r dbe8438d5dca -r 9569e9f35374 include/sampler.h --- a/include/sampler.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/sampler.h Wed Apr 23 10:38:33 2008 +0200 @@ -75,12 +75,12 @@ class DefaultSampler: public Sampler { int phase; - int subsample; + int subsample; // 0,1 = no, 1+ = size of sampling grid int oversample; // 0 = no, 1 = 5x, 2 = 9x, 3 = 16x int sx,sy,osa_samp; // current sample properties public: DefaultSampler(Float *abuffer, int &aw, int &ah): - Sampler(abuffer, aw, ah), phase(-1), subsample(8), oversample(0) {}; + Sampler(abuffer, aw, ah), phase(-1), subsample(4), oversample(0) {}; void init(); int initSampleSet(); bool nextSample(Sample *s); diff -r dbe8438d5dca -r 9569e9f35374 include/scene.h --- a/include/scene.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/scene.h Wed Apr 23 10:38:33 2008 +0200 @@ -30,24 +30,15 @@ #include #include +#include "common.h" #include "sampler.h" #include "noise.h" #include "vector.h" #include "quaternion.h" -/* -triangle intersection alghoritm -chooses are: -TRI_PLUCKER -TRI_BARI -TRI_BARI_PRE -*/ -#if !defined(TRI_PLUCKER) && !defined(TRI_BARI) && !defined(TRI_BARI_PRE) -# define TRI_BARI_PRE -#endif - -using namespace std; - +/** + * ray + */ class Ray { public: @@ -164,167 +155,4 @@ void setSmooth(bool sm) { smooth = sm; }; }; -/** - * shape - */ -class Shape -{ -public: - Material *material; - Shape() {}; - virtual ~Shape() {}; - - // first intersection point - virtual bool intersect(const Ray &ray, Float &dist) const = 0; - - // all intersections (only for CSG) - virtual bool intersect_all(const Ray &ray, Float dist, vector &allts) const = 0; - - // intersection with AABB - virtual bool intersect_bbox(const BBox &bbox) const = 0; - - // normal at point P - virtual const Vector3 normal(const Vector3 &P) const = 0; - - virtual BBox get_bbox() const = 0; -}; - -/** - * list of shapes - */ -class ShapeList: public vector -{ -}; - -/** - * sphere shape - */ -class Sphere: public Shape -{ - Float sqr_radius; - Float inv_radius; -public: - Vector3 center; - Float radius; - - Sphere(const Vector3 &acenter, const Float aradius, Material *amaterial): - sqr_radius(aradius*aradius), inv_radius(1.0f/aradius), - center(acenter), radius(aradius) { material = amaterial; } - bool intersect(const Ray &ray, Float &dist) const; - bool intersect_all(const Ray &ray, Float dist, vector &allts) const; - bool intersect_bbox(const BBox &bbox) const; - const Vector3 normal(const Vector3 &P) const { return (P - center) * inv_radius; }; - BBox get_bbox() const; -}; - -/** - * box shape - */ -class Box: public Shape -{ - Vector3 L; - Vector3 H; -public: - Box(const Vector3 &aL, const Vector3 &aH, Material *amaterial): L(aL), H(aH) - { - for (int i = 0; i < 3; i++) - if (L.cell[i] > H.cell[i]) - swap(L.cell[i], H.cell[i]); - material = amaterial; - }; - bool intersect(const Ray &ray, Float &dist) const; - bool intersect_all(const Ray &ray, Float dist, vector &allts) const { return false; }; - bool intersect_bbox(const BBox &bbox) const; - const Vector3 normal(const Vector3 &P) const; - BBox get_bbox() const { return BBox(L, H); }; -}; - -/** - * triangle vertex - */ -class Vertex -{ -public: - Vector3 P; - Vertex(const Vector3 &aP): P(aP) {}; -}; - -/** - * triangle vertex with normal - */ -class NormalVertex: public Vertex -{ -public: - Vector3 N; - NormalVertex(const NormalVertex *v): Vertex(v->P), N(v->N) {}; - NormalVertex(const Vector3 &aP): Vertex(aP) {}; - NormalVertex(const Vector3 &aP, const Vector3 &aN): Vertex(aP), N(aN) {}; - const Vector3 &getNormal() { return N; }; - void setNormal(const Vector3 &aN) { N = aN; }; -}; - -/** - * triangle shape - */ -class Triangle: public Shape -{ -#ifdef TRI_BARI_PRE - Float nu, nv, nd; - int k; // dominant axis - Float bnu, bnv; - Float cnu, cnv; #endif -#ifdef TRI_BARI - int k; // dominant axis -#endif -#ifdef TRI_PLUCKER - Float pla[6], plb[6], plc[6]; -#endif - Vector3 N; - const Vector3 smooth_normal(const Vector3 &P) const - { -#ifdef TRI_BARI_PRE - const Vector3 &NA = static_cast(A)->N; - const Vector3 &NB = static_cast(B)->N; - const Vector3 &NC = static_cast(C)->N; - static const int modulo3[5] = {0,1,2,0,1}; - register const int ku = modulo3[k+1]; - register const int kv = modulo3[k+2]; - const Float pu = P[ku] - A->P[ku]; - const Float pv = P[kv] - A->P[kv]; - const Float u = pv * bnu + pu * bnv; - const Float v = pu * cnv + pv * cnu; - Vector3 n = NA + u * (NB - NA) + v * (NC - NA); - n.normalize(); - return n; -#else - return N; // not implemented for other algorithms -#endif - }; -public: - Vertex *A, *B, *C; - - Triangle() {}; - Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial); - bool intersect(const Ray &ray, Float &dist) const; - bool intersect_all(const Ray &ray, Float dist, vector &allts) const {return false;}; - bool intersect_bbox(const BBox &bbox) const; - const Vector3 normal(const Vector3 &P) const { return (material->smooth ? smooth_normal(P) : N); }; - const Vector3 getNormal() const { return N; }; - BBox get_bbox() const; -}; - -template class Array -{ - T *array; -public: - Array(int n) { array = new T[n]; }; - ~Array() { delete[] array; }; - const T &operator[](int i) const { return array[i]; }; -}; - -typedef Array VertexArray; -typedef Array NormalVertexArray; -typedef Array TriangleArray; - -#endif diff -r dbe8438d5dca -r 9569e9f35374 include/serialize.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/include/serialize.h Wed Apr 23 10:38:33 2008 +0200 @@ -0,0 +1,57 @@ +/* + * serialize.h: object serialization functions + * + * This file is part of Pyrit Ray Tracer. + * + * Copyright 2008 Radek Brich + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +#ifndef SERIALIZE_H +#define SERIALIZE_H + +#include + +#include "shapes.h" +#include "container.h" +#include "kdtree.h" + +class Indexer +{ + map indexmap; + int index; +public: + Indexer(): indexmap(), index(0) {}; + void reset() { indexmap.clear(); index = 0; }; + bool get(void *o, int &retidx); + const int &operator[](void *o) { return indexmap[o]; }; +}; + +extern Indexer vertex_index, shape_index; + +void resetSerializer(); +Shape *loadShape(istream &st); + +ostream & operator<<(ostream &st, Shape &o); +ostream & operator<<(ostream &st, Vertex &o); +ostream & operator<<(ostream &st, Container &o); +istream & operator>>(istream &st, Vector3 &v); + +#endif diff -r dbe8438d5dca -r 9569e9f35374 include/shapes.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/include/shapes.h Wed Apr 23 10:38:33 2008 +0200 @@ -0,0 +1,217 @@ +/* + * shapes.h: shape classes + * + * This file is part of Pyrit Ray Tracer. + * + * Copyright 2006, 2007, 2008 Radek Brich + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +#ifndef SHAPES_H +#define SHAPES_H + +#include "common.h" +#include "scene.h" + +/* +triangle intersection alghoritm +options are: +TRI_PLUCKER +TRI_BARI +TRI_BARI_PRE +*/ +#if !defined(TRI_PLUCKER) && !defined(TRI_BARI) && !defined(TRI_BARI_PRE) +# define TRI_BARI_PRE +#endif + +/** + * shape + */ +class Shape +{ +public: + Material *material; + Shape() {}; + virtual ~Shape() {}; + + // first intersection point + virtual bool intersect(const Ray &ray, Float &dist) const = 0; + + // all intersections (only for CSG) + virtual bool intersect_all(const Ray &ray, Float dist, vector &allts) const = 0; + + // intersection with AABB + virtual bool intersect_bbox(const BBox &bbox) const = 0; + + // normal at point P + virtual const Vector3 normal(const Vector3 &P) const = 0; + + virtual BBox get_bbox() const = 0; + + virtual ostream & dump(ostream &st) const = 0; +}; + +/** + * list of shapes + */ +typedef vector ShapeList; + +/** + * sphere shape + */ +class Sphere: public Shape +{ + Vector3 center; + Float radius; + + Float sqr_radius; + Float inv_radius; +public: + Sphere(const Vector3 &acenter, const Float aradius, Material *amaterial): + center(acenter), radius(aradius), + sqr_radius(aradius*aradius), inv_radius(1.0f/aradius) + { material = amaterial; } + bool intersect(const Ray &ray, Float &dist) const; + bool intersect_all(const Ray &ray, Float dist, vector &allts) const; + bool intersect_bbox(const BBox &bbox) const; + const Vector3 normal(const Vector3 &P) const { return (P - center) * inv_radius; }; + BBox get_bbox() const; + const Vector3 getCenter() const { return center; }; + const Float getRadius() const { return radius; }; + ostream & dump(ostream &st) const; +}; + +/** + * box shape + */ +class Box: public Shape +{ + Vector3 L; + Vector3 H; +public: + Box(const Vector3 &aL, const Vector3 &aH, Material *amaterial): L(aL), H(aH) + { + for (int i = 0; i < 3; i++) + if (L.cell[i] > H.cell[i]) + swap(L.cell[i], H.cell[i]); + material = amaterial; + }; + bool intersect(const Ray &ray, Float &dist) const; + bool intersect_all(const Ray &ray, Float dist, vector &allts) const { return false; }; + bool intersect_bbox(const BBox &bbox) const; + const Vector3 normal(const Vector3 &P) const; + BBox get_bbox() const { return BBox(L, H); }; + const Vector3 getL() const { return L; }; + const Vector3 getH() const { return H; }; + ostream & dump(ostream &st) const; +}; + +/** + * triangle vertex + */ +class Vertex +{ +public: + Vector3 P; + Vertex(const Vector3 &aP): P(aP) {}; + virtual ostream & dump(ostream &st) const; +}; + +/** + * triangle vertex with normal + */ +class NormalVertex: public Vertex +{ +public: + Vector3 N; + NormalVertex(const NormalVertex *v): Vertex(v->P), N(v->N) {}; + NormalVertex(const Vector3 &aP): Vertex(aP) {}; + NormalVertex(const Vector3 &aP, const Vector3 &aN): Vertex(aP), N(aN) {}; + const Vector3 &getNormal() { return N; }; + void setNormal(const Vector3 &aN) { N = aN; }; + ostream & dump(ostream &st) const; +}; + +/** + * triangle shape + */ +class Triangle: public Shape +{ +#ifdef TRI_BARI_PRE + Float nu, nv, nd; + int k; // dominant axis + Float bnu, bnv; + Float cnu, cnv; +#endif +#ifdef TRI_BARI + int k; // dominant axis +#endif +#ifdef TRI_PLUCKER + Float pla[6], plb[6], plc[6]; +#endif + Vector3 N; + const Vector3 smooth_normal(const Vector3 &P) const + { +#ifdef TRI_BARI_PRE + const Vector3 &NA = static_cast(A)->N; + const Vector3 &NB = static_cast(B)->N; + const Vector3 &NC = static_cast(C)->N; + static const int modulo3[5] = {0,1,2,0,1}; + register const int ku = modulo3[k+1]; + register const int kv = modulo3[k+2]; + const Float pu = P[ku] - A->P[ku]; + const Float pv = P[kv] - A->P[kv]; + const Float u = pv * bnu + pu * bnv; + const Float v = pu * cnv + pv * cnu; + Vector3 n = NA + u * (NB - NA) + v * (NC - NA); + n.normalize(); + return n; +#else + return N; // not implemented for other algorithms +#endif + }; +public: + Vertex *A, *B, *C; + + Triangle() {}; + Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial); + bool intersect(const Ray &ray, Float &dist) const; + bool intersect_all(const Ray &ray, Float dist, vector &allts) const {return false;}; + bool intersect_bbox(const BBox &bbox) const; + const Vector3 normal(const Vector3 &P) const { return (material->smooth ? smooth_normal(P) : N); }; + const Vector3 getNormal() const { return N; }; + BBox get_bbox() const; + ostream & dump(ostream &st) const; +}; + +template class Array +{ + T *array; +public: + Array(int n) { array = new T[n]; }; + ~Array() { delete[] array; }; + const T &operator[](int i) const { return array[i]; }; +}; + +typedef Array VertexArray; +typedef Array NormalVertexArray; +typedef Array TriangleArray; + +#endif diff -r dbe8438d5dca -r 9569e9f35374 include/vector.h --- a/include/vector.h Tue Apr 22 13:33:12 2008 +0200 +++ b/include/vector.h Wed Apr 23 10:38:33 2008 +0200 @@ -30,6 +30,8 @@ #include #include +#include "common.h" + using namespace std; /** @@ -166,10 +168,24 @@ return Vector3(a.x * b.x, a.y * b.y, a.z * b.z); }; - // print + // write friend ostream & operator<<(ostream &st, const Vector3 &v) { - return st << "(" << v.x << ", " << v.y << ", " << v.z << ")"; + return st << "(" << v.x << "," << v.y << "," << v.z << ")"; + }; + + // read + friend istream & operator>>(istream &st, Vector3 &v) + { + char s[10]; + st.getline(s, 10, '('); + st >> v.x; + st.getline(s, 10, ','); + st >> v.y; + st.getline(s, 10, ','); + st >> v.z; + st.getline(s, 10, ')'); + return st; }; }; diff -r dbe8438d5dca -r 9569e9f35374 src/SConscript --- a/src/SConscript Tue Apr 22 13:33:12 2008 +0200 +++ b/src/SConscript Wed Apr 23 10:38:33 2008 +0200 @@ -16,8 +16,9 @@ pyenv.ParseConfig('python-config --includes --libs') sources = [ - 'raytracer.cc', 'scene.cc', 'sampler.cc', - 'container.cc', 'kdtree.cc', 'octree.cc', 'noise.cc'] + 'raytracer.cc', 'scene.cc', 'shapes.cc', 'sampler.cc', + 'container.cc', 'kdtree.cc', 'octree.cc', 'noise.cc', + 'serialize.cc'] objs = [] shared_objs = [] diff -r dbe8438d5dca -r 9569e9f35374 src/container.cc --- a/src/container.cc Tue Apr 22 13:33:12 2008 +0200 +++ b/src/container.cc Wed Apr 23 10:38:33 2008 +0200 @@ -26,11 +26,11 @@ #include "common.h" #include "container.h" +#include "serialize.h" void Container::addShape(Shape* aShape) { - const Float e = 10*Eps; - shapes.push_back(aShape); + const Float e = Eps; if (shapes.size() == 0) { /* initialize bounding box */ bbox = aShape->get_bbox(); @@ -46,6 +46,7 @@ if (shapebb.H.y + e > bbox.H.y) bbox.H.y = shapebb.H.y + e; if (shapebb.H.z + e > bbox.H.z) bbox.H.z = shapebb.H.z + e; } + shapes.push_back(aShape); }; Shape *Container::nearest_intersection(const Shape *origin_shape, const Ray &ray, @@ -58,3 +59,16 @@ nearest_shape = *shape; return nearest_shape; } + +ostream & Container::dump(ostream &st) +{ + st << "(container," << shapes.size(); + ShapeList::iterator shape; + for (shape = shapes.begin(); shape != shapes.end(); shape++) + { + int idx; + if (!shape_index.get(*shape, idx)) + st << "," << **shape; + } + return st << ")"; +} diff -r dbe8438d5dca -r 9569e9f35374 src/kdtree.cc --- a/src/kdtree.cc Tue Apr 22 13:33:12 2008 +0200 +++ b/src/kdtree.cc Wed Apr 23 10:38:33 2008 +0200 @@ -26,9 +26,11 @@ #include #include +#include +#include -#include "common.h" #include "kdtree.h" +#include "serialize.h" class ShapeBound { @@ -265,7 +267,7 @@ } else { /* (stack[entry].pb[axis] > splitVal) */ - if (splitVal < stack[exit].pb[axis]) + if (stack[exit].pb[axis] > splitVal) { /* case P1, P2, P3, and N5 */ node = node->getRightChild(); @@ -305,7 +307,7 @@ ShapeList::iterator shape; for (shape = node->getShapes()->begin(); shape != node->getShapes()->end(); shape++) if (*shape != origin_shape && (*shape)->intersect(ray, dist) - && dist >= stack[entry].t) + && dist >= stack[entry].t - Eps) { nearest_shape = *shape; nearest_distance = dist; @@ -326,27 +328,118 @@ return NULL; } -// this should save whole kd-tree with triangles distributed into leaves -void KdTree::save(ostream &str, KdNode *node) +ostream & operator<<(ostream &st, KdNode &node) +{ + if (node.isLeaf()) + { + st << "(l," << node.getShapes()->size(); + ShapeList::iterator shape; + for (shape = node.getShapes()->begin(); shape != node.getShapes()->end(); shape++) + st << "," << shape_index[*shape]; + st << ")"; + } + else + { + st << "(s," << (char)('x'+node.getAxis()) << "," << node.getSplit() << ","; + st << *node.getLeftChild() << ","; + st << *node.getRightChild() << ")"; + } + return st; +} + +ostream & KdTree::dump(ostream &st) { if (!built) - return; - if (node == NULL) - node = root; - if (node->isLeaf()) - str << "(leaf: " << node->getShapes()->size() << " shapes)"; - else + return Container::dump(st); + + st << "(kdtree," << shapes.size(); + ShapeList::iterator shape; + for (shape = shapes.begin(); shape != shapes.end(); shape++) + { + int idx; + if (!shape_index.get(*shape, idx)) + st << "," << endl << **shape; + } + return st << "," << endl << *getRootNode() << ")"; +} + +void KdTree::recursive_load(istream &st, KdNode *node) +{ + string s; + istringstream is; + + getline(st, s, ','); + trim(s); + + if (s.compare("(s") == 0) { - str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L="; - save(str, node->getLeftChild()); - str << "; R="; - save(str, node->getRightChild()); - str << ";)"; + // split + int axis; + Float split; + + delete node->getShapes(); + node->setChildren(new KdNode[2]); + + getline(st, s, ','); + axis = s.c_str()[0]-'x'; + node->setAxis(axis); + + st >> split; + getline(st, s, ','); + node->setSplit(split); + + recursive_load(st, node->getLeftChild()); + getline(st, s, ','); + recursive_load(st, node->getRightChild()); + getline(st, s, ')'); + } + + if (s.compare("(l") == 0) + { + // leaf + int count, idx; + + node->setLeaf(); + + st >> count; + for (int i = 0; i < count; i++) + { + getline(st, s, ','); + st >> idx; + node->addShape(shapes[idx]); + } + getline(st, s, ')'); } } -// load kd-tree from file/stream -void KdTree::load(istream &str, KdNode *node) +istream & KdTree::load(istream &st) { + string s; + istringstream is; + getline(st, s, ','); + if (s.compare("(kdtree") != 0) + return st; + + shapes.clear(); + if (root) delete root; + root = new KdNode(); + + getline(st, s, ','); + int shape_count; + is.str(s); + is >> shape_count; + + Shape *shape; + for (int i = 0; i < shape_count; i++) + { + shape = loadShape(st); + Container::addShape(shape); + getline(st, s, ','); + } + + recursive_load(st, root); + + built = true; + return st; } diff -r dbe8438d5dca -r 9569e9f35374 src/raytracermodule.cc --- a/src/raytracermodule.cc Tue Apr 22 13:33:12 2008 +0200 +++ b/src/raytracermodule.cc Wed Apr 23 10:38:33 2008 +0200 @@ -28,7 +28,6 @@ #include #include "raytracer.h" -#include "scene.h" #include "octree.h" //=========================== Light Source Object =========================== diff -r dbe8438d5dca -r 9569e9f35374 src/scene.cc --- a/src/scene.cc Tue Apr 22 13:33:12 2008 +0200 +++ b/src/scene.cc Wed Apr 23 10:38:33 2008 +0200 @@ -24,9 +24,6 @@ * THE SOFTWARE. */ -#include - -#include "common.h" #include "scene.h" void Camera::rotate(const Quaternion &q) @@ -110,474 +107,3 @@ b = tfar; return true; } - -bool Sphere::intersect(const Ray &ray, Float &dist) const -{ - Vector3 V = ray.o - center; - register Float d = -dot(V, ray.dir); - register Float Det = d * d - (dot(V,V) - sqr_radius); - register Float t1,t2; - if (Det > 0) { - Det = sqrtf(Det); - t1 = d - Det; - t2 = d + Det; - if (t1 > 0) - { - if (t1 < dist) - { - dist = t1; - return true; - } - } - else if (t2 > 0 && t2 < dist) - { - dist = t2; - return true; - } - } - return false; -} - -/* if there should be CSG sometimes, this may be needed... */ -bool Sphere::intersect_all(const Ray &ray, Float dist, vector &allts) const -{ - //allts = new vector(); - - Vector3 V = ((Ray)ray).o - center; - Float Vd = - dot(V, ray.dir); - Float Det = Vd * Vd - (dot(V,V) - sqr_radius); - - if (Det > 0) { - Det = sqrtf(Det); - Float t1 = Vd - Det; - Float t2 = Vd + Det; - if (t1 < 0) - { - if (t2 > 0) - { - // ray from inside of the sphere - allts.push_back(0.0); - allts.push_back(t2); - return true; - } - else - return false; - } - else - { - allts.push_back(t1); - allts.push_back(t2); - return true; - } - } - return false; -} - -bool Sphere::intersect_bbox(const BBox &bbox) const -{ - register float dmin = 0; - for (int i = 0; i < 3; i++) - { - if (center[i] < bbox.L[i]) - dmin += (center[i] - bbox.L[i])*(center[i] - bbox.L[i]); - else - if (center[i] > bbox.H[i]) - dmin += (center[i] - bbox.H[i])*(center[i] - bbox.H[i]); - } - if (dmin <= sqr_radius) - return true; - return false; -}; - -BBox Sphere::get_bbox() const -{ - return BBox(center - radius, center + radius); -} - -bool Box::intersect(const Ray &ray, Float &dist) const -{ - register Float tnear = -Inf; - register Float tfar = Inf; - register Float t1, t2; - - for (int i = 0; i < 3; i++) - { - if (ray.dir[i] == 0) { - /* ray is parallel to these planes */ - if (ray.o[i] < L[i] || ray.o[i] > H[i]) - return false; - } - else - { - /* compute the intersection distance of the planes */ - t1 = (L[i] - ray.o[i]) / ray.dir[i]; - t2 = (H[i] - ray.o[i]) / ray.dir[i]; - - if (t1 > t2) - swap(t1, t2); - - if (t1 > tnear) - tnear = t1; /* want largest Tnear */ - if (t2 < tfar) - tfar = t2; /* want smallest Tfar */ - if (tnear > tfar || tfar < 0) - return false; /* box missed; box is behind ray */ - } - } - - if (tnear < dist) - { - dist = tnear; - return true; - } - return false; -} - -bool Box::intersect_bbox(const BBox &bbox) const -{ - return ( - H.x > bbox.L.x && L.x < bbox.H.x && - H.y > bbox.L.y && L.y < bbox.H.y && - H.z > bbox.L.z && L.z < bbox.H.z); -} - -const Vector3 Box::normal(const Vector3 &P) const -{ - register Vector3 l = P - L; - register Vector3 h = H - P; - - if (l.x < h.x) - h.x = -1; - else - { - l.x = h.x; - h.x = +1; - } - - if (l.y < h.y) - h.y = -1; - else - { - l.y = h.y; - h.y = +1; - } - - if (l.z < h.z) - h.z = -1; - else - { - l.z = h.z; - h.z = +1; - } - - if (l.x > l.y) - { - h.x = 0; - if (l.y > l.z) - h.y = 0; - else - h.z = 0; - } - else - { - h.y = 0; - if (l.x > l.z) - h.x = 0; - else - h.z = 0; - } - return h; -} - -#ifdef TRI_PLUCKER -inline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl) -{ - pl[0] = p.x*q.y - q.x*p.y; - pl[1] = p.x*q.z - q.x*p.z; - pl[2] = p.x - q.x; - pl[3] = p.y*q.z - q.y*p.z; - pl[4] = p.z - q.z; - pl[5] = q.y - p.y; -} - -inline Float Side(const Float* pla, const Float* plb) -{ - return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2]; -} -#endif - -Triangle::Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial) - : A(aA), B(aB), C(aC) -{ - material = amaterial; - - const Vector3 c = B->P - A->P; - const Vector3 b = C->P - A->P; - - N = cross(c, b); - N.normalize(); - -#ifdef TRI_PLUCKER - Plucker(B->P,C->P,pla); - Plucker(C->P,A->P,plb); - Plucker(A->P,B->P,plc); -#endif - -#if defined(TRI_BARI) || defined(TRI_BARI_PRE) - if (fabsf(N.x) > fabsf(N.y)) - { - if (fabsf(N.x) > fabsf(N.z)) - k = 0; - else - k = 2; - } - else - { - if (fabsf(N.y) > fabsf(N.z)) - k = 1; - else - k = 2; - } -#endif -#ifdef TRI_BARI_PRE - int u = (k + 1) % 3; - int v = (k + 2) % 3; - - Float krec = 1.0 / N[k]; - nu = N[u] * krec; - nv = N[v] * krec; - nd = dot(N, A->P) * krec; - - // first line equation - Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]); - bnu = b[u] * reci; - bnv = -b[v] * reci; - - // second line equation - cnu = -c[u] * reci; - cnv = c[v] * reci; -#endif -} - -bool Triangle::intersect(const Ray &ray, Float &dist) const -{ -#ifdef TRI_PLUCKER - Float plr[6]; - Plucker(ray.o, ray.o+ray.dir, plr); - const bool side0 = Side(plr, pla) >= 0.0; - const bool side1 = Side(plr, plb) >= 0.0; - if (side0 != side1) - return false; - const bool side2 = Side(plr, plc) >= 0.0; - if (side0 != side2) - return false; - const Float t = - dot( (ray.o - A->P), N) / dot(ray.dir,N); - if(t <= Eps || t >= dist) - return false; - dist = t; - return true; -#endif - -#if defined(TRI_BARI) || defined(TRI_BARI_PRE) - static const int modulo3[5] = {0,1,2,0,1}; - const Vector3 &O = ray.o; - const Vector3 &D = ray.dir; - register const int u = modulo3[k+1]; - register const int v = modulo3[k+2]; -#endif -#ifdef TRI_BARI_PRE - const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]); - - if (t >= dist || t < Eps) - return false; - - const Float hu = O[u] + t * D[u] - A->P[u]; - const Float hv = O[v] + t * D[v] - A->P[v]; - const Float beta = hv * bnu + hu * bnv; - - if (beta < 0.) - return false; - - const Float gamma = hu * cnv + hv * cnu; - if (gamma < 0. || beta + gamma > 1.) - return false; - - dist = t; - return true; -#endif - -#ifdef TRI_BARI - // original barycentric coordinates based intesection - // not optimized, just for reference - const Vector3 c = B - A; - const Vector3 b = C - A; - // distance test - const Float t = - dot( (O-A), N) / dot(D,N); - if (t < Eps || t > dist) - return false; - - // calc hitpoint - const Float Hu = O[u] + t * D[u] - A[u]; - const Float Hv = O[v] + t * D[v] - A[v]; - const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]); - if (beta < 0) - return false; - const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]); - if (gamma < 0) - return false; - if (beta+gamma > 1) - return false; - dist = t; - return true; -#endif -} - -bool Triangle::intersect_bbox(const BBox &bbox) const -{ - const Vector3 boxcenter = (bbox.L+bbox.H)*0.5; - const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5; - const Vector3 v0 = A->P - boxcenter; - const Vector3 v1 = B->P - boxcenter; - const Vector3 v2 = C->P - boxcenter; - const Vector3 e0 = v1-v0; - const Vector3 e1 = v2-v1; - const Vector3 e2 = v0-v2; - - Float fex = fabsf(e0.x); - Float fey = fabsf(e0.y); - Float fez = fabsf(e0.z); - - Float p0,p1,p2,min,max,rad; - - p0 = e0.z*v0.y - e0.y*v0.z; - p2 = e0.z*v2.y - e0.y*v2.z; - if(p0rad || max<-rad) return false; - - p0 = -e0.z*v0.x + e0.x*v0.z; - p2 = -e0.z*v2.x + e0.x*v2.z; - if(p0rad || max<-rad) return false; - - p1 = e0.y*v1.x - e0.x*v1.y; - p2 = e0.y*v2.x - e0.x*v2.y; - if(p2rad || max<-rad) return false; - - fex = fabsf(e1.x); - fey = fabsf(e1.y); - fez = fabsf(e1.z); - - p0 = e1.z*v0.y - e1.y*v0.z; - p2 = e1.z*v2.y - e1.y*v2.z; - if(p0rad || max<-rad) return false; - - p0 = -e1.z*v0.x + e1.x*v0.z; - p2 = -e1.z*v2.x + e1.x*v2.z; - if(p0rad || max<-rad) return false; - - p0 = e1.y*v0.x - e1.x*v0.y; - p1 = e1.y*v1.x - e1.x*v1.y; - if(p0rad || max<-rad) return false; - - fex = fabsf(e2.x); - fey = fabsf(e2.y); - fez = fabsf(e2.z); - - p0 = e2.z*v0.y - e2.y*v0.z; - p1 = e2.z*v1.y - e2.y*v1.z; - if(p0rad || max<-rad) return false; - - p0 = -e2.z*v0.x + e2.x*v0.z; - p1 = -e2.z*v1.x + e2.x*v1.z; - if(p0rad || max<-rad) return false; - - p1 = e2.y*v1.x - e2.x*v1.y; - p2 = e2.y*v2.x - e2.x*v2.y; - if(p2rad || max<-rad) return false; - - /* test overlap in the {x,y,z}-directions */ - /* test in X-direction */ - min = v0.x; - if (v1.x < min) min = v1.x; - if (v2.x < min) min = v2.x; - max = v0.x; - if (v1.x > max) max = v1.x; - if (v2.x > max) max = v2.x; - if(min>boxhalfsize.x || max<-boxhalfsize.x) return false; - - /* test in Y-direction */ - min = v0.y; - if (v1.y < min) min = v1.y; - if (v2.y < min) min = v2.y; - max = v0.y; - if (v1.y > max) max = v1.y; - if (v2.y > max) max = v2.y; - if(min>boxhalfsize.y || max<-boxhalfsize.y) return false; - - /* test in Z-direction */ - min = v0.z; - if (v1.z < min) min = v1.z; - if (v2.z < min) min = v2.z; - max = v0.z; - if (v1.z > max) max = v1.z; - if (v2.z > max) max = v2.z; - if(min>boxhalfsize.z || max<-boxhalfsize.z) return false; - - /* test if the box intersects the plane of the triangle */ - Vector3 vmin,vmax; - Float v; - for(int q=0;q<3;q++) - { - v=v0[q]; - if(N[q]>0.0f) - { - vmin.cell[q]=-boxhalfsize[q] - v; - vmax.cell[q]= boxhalfsize[q] - v; - } - else - { - vmin.cell[q]= boxhalfsize[q] - v; - vmax.cell[q]=-boxhalfsize[q] - v; - } - } - if(dot(N,vmin)>0.0f) return false; - if(dot(N,vmax)>=0.0f) return true; - - return false; -} - -BBox Triangle::get_bbox() const -{ - BBox bbox = BBox(); - bbox.L = A->P; - if (B->P.x < bbox.L.x) bbox.L.x = B->P.x; - if (C->P.x < bbox.L.x) bbox.L.x = C->P.x; - if (B->P.y < bbox.L.y) bbox.L.y = B->P.y; - if (C->P.y < bbox.L.y) bbox.L.y = C->P.y; - if (B->P.z < bbox.L.z) bbox.L.z = B->P.z; - if (C->P.z < bbox.L.z) bbox.L.z = C->P.z; - bbox.H = A->P; - if (B->P.x > bbox.H.x) bbox.H.x = B->P.x; - if (C->P.x > bbox.H.x) bbox.H.x = C->P.x; - if (B->P.y > bbox.H.y) bbox.H.y = B->P.y; - if (C->P.y > bbox.H.y) bbox.H.y = C->P.y; - if (B->P.z > bbox.H.z) bbox.H.z = B->P.z; - if (C->P.z > bbox.H.z) bbox.H.z = C->P.z; - return bbox; -}; diff -r dbe8438d5dca -r 9569e9f35374 src/serialize.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/serialize.cc Wed Apr 23 10:38:33 2008 +0200 @@ -0,0 +1,97 @@ +/* + * serialize.cc: object serialization functions + * + * This file is part of Pyrit Ray Tracer. + * + * Copyright 2008 Radek Brich + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +#include "serialize.h" +#include +#include + +Indexer vertex_index, shape_index; + +void resetSerializer() +{ + vertex_index.reset(); + shape_index.reset(); +} + +bool Indexer::get(void *o, int &retidx) +{ + map ::iterator i; + i = indexmap.find(o); + if (i == indexmap.end()) + { + retidx = index++; + indexmap[o] = retidx; + return false; + } + else + { + retidx = i->second; + return true; + } +} + +Shape *loadShape(istream &st) +{ + string s; + istringstream is; + getline(st, s, ','); + trim(s); + if (s.compare("(box") == 0) + { + Vector3 L,H; + st >> L; + getline(st, s, ','); + st >> H; + getline(st, s, ')'); + return new Box(L, H, new Material(Colour(1,1,1))); + } + if (s.compare("(sphere") == 0) + { + Vector3 center; + Float radius; + st >> center; + getline(st, s, ','); + st >> radius; + getline(st, s, ')'); + return new Sphere(center, radius, new Material(Colour(1,1,1))); + } + return NULL; +} + +ostream & operator<<(ostream &st, Shape &o) +{ + return o.dump(st); +} + +ostream & operator<<(ostream &st, Vertex &o) +{ + return o.dump(st); +} + +ostream & operator<<(ostream &st, Container &o) +{ + return o.dump(st); +} diff -r dbe8438d5dca -r 9569e9f35374 src/shapes.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/shapes.cc Wed Apr 23 10:38:33 2008 +0200 @@ -0,0 +1,531 @@ +/* + * shapes.cc: shape classes + * + * This file is part of Pyrit Ray Tracer. + * + * Copyright 2006, 2007, 2008 Radek Brich + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +#include "shapes.h" +#include "serialize.h" + +bool Sphere::intersect(const Ray &ray, Float &dist) const +{ + Vector3 V = ray.o - center; + register Float d = -dot(V, ray.dir); + register Float Det = d * d - (dot(V,V) - sqr_radius); + register Float t1,t2; + if (Det > 0) { + Det = sqrtf(Det); + t1 = d - Det; + t2 = d + Det; + if (t1 > 0) + { + if (t1 < dist) + { + dist = t1; + return true; + } + } + else if (t2 > 0 && t2 < dist) + { + dist = t2; + return true; + } + } + return false; +} + +/* if there should be CSG sometimes, this may be needed... */ +bool Sphere::intersect_all(const Ray &ray, Float dist, vector &allts) const +{ + //allts = new vector(); + + Vector3 V = ((Ray)ray).o - center; + Float Vd = - dot(V, ray.dir); + Float Det = Vd * Vd - (dot(V,V) - sqr_radius); + + if (Det > 0) { + Det = sqrtf(Det); + Float t1 = Vd - Det; + Float t2 = Vd + Det; + if (t1 < 0) + { + if (t2 > 0) + { + // ray from inside of the sphere + allts.push_back(0.0); + allts.push_back(t2); + return true; + } + else + return false; + } + else + { + allts.push_back(t1); + allts.push_back(t2); + return true; + } + } + return false; +} + +bool Sphere::intersect_bbox(const BBox &bbox) const +{ + register float dmin = 0; + for (int i = 0; i < 3; i++) + { + if (center[i] < bbox.L[i]) + dmin += (center[i] - bbox.L[i])*(center[i] - bbox.L[i]); + else + if (center[i] > bbox.H[i]) + dmin += (center[i] - bbox.H[i])*(center[i] - bbox.H[i]); + } + if (dmin <= sqr_radius) + return true; + return false; +}; + +BBox Sphere::get_bbox() const +{ + return BBox(center - radius, center + radius); +} + +bool Box::intersect(const Ray &ray, Float &dist) const +{ + register Float tnear = -Inf; + register Float tfar = Inf; + register Float t1, t2; + + for (int i = 0; i < 3; i++) + { + if (ray.dir[i] == 0) { + /* ray is parallel to these planes */ + if (ray.o[i] < L[i] || ray.o[i] > H[i]) + return false; + } + else + { + /* compute the intersection distance of the planes */ + t1 = (L[i] - ray.o[i]) / ray.dir[i]; + t2 = (H[i] - ray.o[i]) / ray.dir[i]; + + if (t1 > t2) + swap(t1, t2); + + if (t1 > tnear) + tnear = t1; /* want largest Tnear */ + if (t2 < tfar) + tfar = t2; /* want smallest Tfar */ + if (tnear > tfar || tfar < 0) + return false; /* box missed; box is behind ray */ + } + } + + if (tnear < dist) + { + dist = tnear; + return true; + } + return false; +} + +bool Box::intersect_bbox(const BBox &bbox) const +{ + return ( + H.x > bbox.L.x && L.x < bbox.H.x && + H.y > bbox.L.y && L.y < bbox.H.y && + H.z > bbox.L.z && L.z < bbox.H.z); +} + +const Vector3 Box::normal(const Vector3 &P) const +{ + register Vector3 l = P - L; + register Vector3 h = H - P; + + if (l.x < h.x) + h.x = -1; + else + { + l.x = h.x; + h.x = +1; + } + + if (l.y < h.y) + h.y = -1; + else + { + l.y = h.y; + h.y = +1; + } + + if (l.z < h.z) + h.z = -1; + else + { + l.z = h.z; + h.z = +1; + } + + if (l.x > l.y) + { + h.x = 0; + if (l.y > l.z) + h.y = 0; + else + h.z = 0; + } + else + { + h.y = 0; + if (l.x > l.z) + h.x = 0; + else + h.z = 0; + } + return h; +} + +#ifdef TRI_PLUCKER +inline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl) +{ + pl[0] = p.x*q.y - q.x*p.y; + pl[1] = p.x*q.z - q.x*p.z; + pl[2] = p.x - q.x; + pl[3] = p.y*q.z - q.y*p.z; + pl[4] = p.z - q.z; + pl[5] = q.y - p.y; +} + +inline Float Side(const Float* pla, const Float* plb) +{ + return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2]; +} +#endif + +Triangle::Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial) + : A(aA), B(aB), C(aC) +{ + material = amaterial; + + const Vector3 c = B->P - A->P; + const Vector3 b = C->P - A->P; + + N = cross(c, b); + N.normalize(); + +#ifdef TRI_PLUCKER + Plucker(B->P,C->P,pla); + Plucker(C->P,A->P,plb); + Plucker(A->P,B->P,plc); +#endif + +#if defined(TRI_BARI) || defined(TRI_BARI_PRE) + if (fabsf(N.x) > fabsf(N.y)) + { + if (fabsf(N.x) > fabsf(N.z)) + k = 0; + else + k = 2; + } + else + { + if (fabsf(N.y) > fabsf(N.z)) + k = 1; + else + k = 2; + } +#endif +#ifdef TRI_BARI_PRE + int u = (k + 1) % 3; + int v = (k + 2) % 3; + + Float krec = 1.0 / N[k]; + nu = N[u] * krec; + nv = N[v] * krec; + nd = dot(N, A->P) * krec; + + // first line equation + Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]); + bnu = b[u] * reci; + bnv = -b[v] * reci; + + // second line equation + cnu = -c[u] * reci; + cnv = c[v] * reci; +#endif +} + +bool Triangle::intersect(const Ray &ray, Float &dist) const +{ +#ifdef TRI_PLUCKER + Float plr[6]; + Plucker(ray.o, ray.o+ray.dir, plr); + const bool side0 = Side(plr, pla) >= 0.0; + const bool side1 = Side(plr, plb) >= 0.0; + if (side0 != side1) + return false; + const bool side2 = Side(plr, plc) >= 0.0; + if (side0 != side2) + return false; + const Float t = - dot( (ray.o - A->P), N) / dot(ray.dir,N); + if(t <= Eps || t >= dist) + return false; + dist = t; + return true; +#endif + +#if defined(TRI_BARI) || defined(TRI_BARI_PRE) + static const int modulo3[5] = {0,1,2,0,1}; + const Vector3 &O = ray.o; + const Vector3 &D = ray.dir; + register const int u = modulo3[k+1]; + register const int v = modulo3[k+2]; +#endif +#ifdef TRI_BARI_PRE + const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]); + + if (t >= dist || t < Eps) + return false; + + const Float hu = O[u] + t * D[u] - A->P[u]; + const Float hv = O[v] + t * D[v] - A->P[v]; + const Float beta = hv * bnu + hu * bnv; + + if (beta < 0.) + return false; + + const Float gamma = hu * cnv + hv * cnu; + if (gamma < 0. || beta + gamma > 1.) + return false; + + dist = t; + return true; +#endif + +#ifdef TRI_BARI + // original barycentric coordinates based intesection + // not optimized, just for reference + const Vector3 c = B - A; + const Vector3 b = C - A; + // distance test + const Float t = - dot( (O-A), N) / dot(D,N); + if (t < Eps || t > dist) + return false; + + // calc hitpoint + const Float Hu = O[u] + t * D[u] - A[u]; + const Float Hv = O[v] + t * D[v] - A[v]; + const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]); + if (beta < 0) + return false; + const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]); + if (gamma < 0) + return false; + if (beta+gamma > 1) + return false; + dist = t; + return true; +#endif +} + +bool Triangle::intersect_bbox(const BBox &bbox) const +{ + const Vector3 boxcenter = (bbox.L+bbox.H)*0.5; + const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5; + const Vector3 v0 = A->P - boxcenter; + const Vector3 v1 = B->P - boxcenter; + const Vector3 v2 = C->P - boxcenter; + const Vector3 e0 = v1-v0; + const Vector3 e1 = v2-v1; + const Vector3 e2 = v0-v2; + + Float fex = fabsf(e0.x); + Float fey = fabsf(e0.y); + Float fez = fabsf(e0.z); + + Float p0,p1,p2,min,max,rad; + + p0 = e0.z*v0.y - e0.y*v0.z; + p2 = e0.z*v2.y - e0.y*v2.z; + if(p0rad || max<-rad) return false; + + p0 = -e0.z*v0.x + e0.x*v0.z; + p2 = -e0.z*v2.x + e0.x*v2.z; + if(p0rad || max<-rad) return false; + + p1 = e0.y*v1.x - e0.x*v1.y; + p2 = e0.y*v2.x - e0.x*v2.y; + if(p2rad || max<-rad) return false; + + fex = fabsf(e1.x); + fey = fabsf(e1.y); + fez = fabsf(e1.z); + + p0 = e1.z*v0.y - e1.y*v0.z; + p2 = e1.z*v2.y - e1.y*v2.z; + if(p0rad || max<-rad) return false; + + p0 = -e1.z*v0.x + e1.x*v0.z; + p2 = -e1.z*v2.x + e1.x*v2.z; + if(p0rad || max<-rad) return false; + + p0 = e1.y*v0.x - e1.x*v0.y; + p1 = e1.y*v1.x - e1.x*v1.y; + if(p0rad || max<-rad) return false; + + fex = fabsf(e2.x); + fey = fabsf(e2.y); + fez = fabsf(e2.z); + + p0 = e2.z*v0.y - e2.y*v0.z; + p1 = e2.z*v1.y - e2.y*v1.z; + if(p0rad || max<-rad) return false; + + p0 = -e2.z*v0.x + e2.x*v0.z; + p1 = -e2.z*v1.x + e2.x*v1.z; + if(p0rad || max<-rad) return false; + + p1 = e2.y*v1.x - e2.x*v1.y; + p2 = e2.y*v2.x - e2.x*v2.y; + if(p2rad || max<-rad) return false; + + /* test overlap in the {x,y,z}-directions */ + /* test in X-direction */ + min = v0.x; + if (v1.x < min) min = v1.x; + if (v2.x < min) min = v2.x; + max = v0.x; + if (v1.x > max) max = v1.x; + if (v2.x > max) max = v2.x; + if(min>boxhalfsize.x || max<-boxhalfsize.x) return false; + + /* test in Y-direction */ + min = v0.y; + if (v1.y < min) min = v1.y; + if (v2.y < min) min = v2.y; + max = v0.y; + if (v1.y > max) max = v1.y; + if (v2.y > max) max = v2.y; + if(min>boxhalfsize.y || max<-boxhalfsize.y) return false; + + /* test in Z-direction */ + min = v0.z; + if (v1.z < min) min = v1.z; + if (v2.z < min) min = v2.z; + max = v0.z; + if (v1.z > max) max = v1.z; + if (v2.z > max) max = v2.z; + if(min>boxhalfsize.z || max<-boxhalfsize.z) return false; + + /* test if the box intersects the plane of the triangle */ + Vector3 vmin,vmax; + Float v; + for(int q=0;q<3;q++) + { + v=v0[q]; + if(N[q]>0.0f) + { + vmin.cell[q]=-boxhalfsize[q] - v; + vmax.cell[q]= boxhalfsize[q] - v; + } + else + { + vmin.cell[q]= boxhalfsize[q] - v; + vmax.cell[q]=-boxhalfsize[q] - v; + } + } + if(dot(N,vmin)>0.0f) return false; + if(dot(N,vmax)>=0.0f) return true; + + return false; +} + +BBox Triangle::get_bbox() const +{ + BBox bbox = BBox(); + bbox.L = A->P; + if (B->P.x < bbox.L.x) bbox.L.x = B->P.x; + if (C->P.x < bbox.L.x) bbox.L.x = C->P.x; + if (B->P.y < bbox.L.y) bbox.L.y = B->P.y; + if (C->P.y < bbox.L.y) bbox.L.y = C->P.y; + if (B->P.z < bbox.L.z) bbox.L.z = B->P.z; + if (C->P.z < bbox.L.z) bbox.L.z = C->P.z; + bbox.H = A->P; + if (B->P.x > bbox.H.x) bbox.H.x = B->P.x; + if (C->P.x > bbox.H.x) bbox.H.x = C->P.x; + if (B->P.y > bbox.H.y) bbox.H.y = B->P.y; + if (C->P.y > bbox.H.y) bbox.H.y = C->P.y; + if (B->P.z > bbox.H.z) bbox.H.z = B->P.z; + if (C->P.z > bbox.H.z) bbox.H.z = C->P.z; + return bbox; +} + +ostream & Sphere::dump(ostream &st) const +{ + return st << "(sphere," << center << "," << radius << ")"; +} + +ostream & Box::dump(ostream &st) const +{ + return st << "(box," << L << "," << H << ")"; +} + +ostream & Vertex::dump(ostream &st) const +{ + return st << "(v," << P << ")"; +} + +ostream & NormalVertex::dump(ostream &st) const +{ + return st << "(vn," << P << "," << N << ")"; +} + +ostream & Triangle::dump(ostream &st) const +{ + int idxA, idxB, idxC; + if (!vertex_index.get(A, idxA)) + st << *A << "," << endl; + if (!vertex_index.get(B, idxB)) + st << *B << "," << endl; + if (!vertex_index.get(C, idxC)) + st << *C << "," << endl; + return st << "(t," << idxA << "," << idxB << "," << idxC << ")"; +}