gKit2 light
bvh.cpp
1 
2 #include <vector>
3 #include <algorithm>
4 #include <cassert>
5 
6 #include "bvh.h"
7 
8 
11 bool intersect( const Triangle *triangle, const Ray& ray, const float htmax, float &t, float &tu, float &tv )
12 {
13  /* begin calculating determinant - also used to calculate U parameter */
14  Vector ac= Vector(triangle->a, triangle->c);
15  const Vector pvec= cross(ray.direction, ac);
16 
17  /* if determinant is near zero, ray lies in plane of triangle */
18  Vector ab= Vector(triangle->a, triangle->b);
19  const float det= dot(ab, pvec);
20  if (det > -0.0001f && det < 0.0001f) return false;
21 
22  const float inv_det= 1.0f / det;
23 
24  /* calculate distance from vert0 to ray origin */
25  const Vector tvec= Vector(triangle->a, ray.origin);
26 
27  /* calculate U parameter and test bounds */
28  const float u= dot(tvec, pvec) * inv_det;
29  if(u < 0.0f || u > 1.0f) return false;
30 
31  /* prepare to test V parameter */
32  const Vector qvec= cross(tvec, ab);
33 
34  /* calculate V parameter and test bounds */
35  const float v= dot(ray.direction, qvec) * inv_det;
36  if(v < 0.0f || u + v > 1.0f) return false;
37 
38  /* calculate t, ray intersects triangle */
39  float hit= dot(ac, qvec) * inv_det;
40 
41  // ne renvoie vrai que si l'intersection est valide (comprise entre tmin / 0 et tmax du rayon)
42  t= hit;
43  tu= u;
44  tv= v;
45  return (hit > 0 && hit < htmax);
46 }
47 
66 bool intersect( const BBox& box, const Ray &ray, const float htmax )
67 {
68  // xslab
69  float tmin= (box.min.x - ray.origin.x) / ray.direction.x;
70  float tmax= (box.max.x - ray.origin.x) / ray.direction.x;
71  if(tmax < tmin) std::swap(tmin, tmax);
72 
73  // y slab
74  float tymin= (box.min.y - ray.origin.y) / ray.direction.y;
75  float tymax= (box.max.y - ray.origin.y) / ray.direction.y;
76  if(tymax < tymin) std::swap(tymin, tymax);
77 
78  if((tmin > tymax) || (tymin > tmax)) return false;
79  if(tymin > tmin) tmin= tymin;
80  if(tymax < tmax) tmax= tymax;
81 
82  // z slab
83  float tzmin= (box.min.z - ray.origin.z) / ray.direction.z;
84  float tzmax= (box.max.z - ray.origin.z) / ray.direction.z;
85  if(tzmax < tzmin) std::swap(tzmin, tzmax);
86 
87  if((tmin > tzmax) || (tzmin > tmax)) return false;
88  if(tzmin > tmin) tmin= tzmin;
89  if(tzmax < tmax) tmax= tzmax;
90 
91  // ne renvoie vrai que si l'intersection est valide
92  return (tmin < tmax && tmax > 0 && tmin < htmax);
93 }
94 
95 
96 BVHNode *create_node( const BBox& box, BVHNode *left, BVHNode *right )
97 {
98  BVHNode *node= new BVHNode;
99  node->left= left;
100  node->right= right;
101  node->triangle= nullptr;
102  node->box= box;
103  return node;
104 }
105 
106 BVHNode *create_leaf( const BBox& box, Triangle *triangle )
107 {
108  BVHNode *leaf= new BVHNode;
109  leaf->left= nullptr;
110  leaf->right= nullptr;
111  leaf->triangle= triangle;
112  leaf->box= box;
113  return leaf;
114 }
115 
116 
117 struct CutLess
118 {
119  CutLess( const int a, const float c ) : axis(a), cut(c) {}
120 
121  int axis;
122  float cut;
123 
124  bool operator() ( const Triangle& t )
125  {
126  // construit la bbox du triangle
127  BBox box;
128  bbox_insert(box, t.a);
129  bbox_insert(box, t.b);
130  bbox_insert(box, t.c);
131 
132  if(axis == 0 && box.min.x < cut) return true;
133  else if(axis == 1 && box.min.y < cut) return true;
134  else if(axis == 2 && box.min.z < cut) return true;
135  return false;
136  }
137 };
138 
139 BVHNode *build_node( std::vector<Triangle>& triangles, const unsigned int begin, const unsigned int end )
140 {
141  if(begin == end)
142  return nullptr; // plus de triangles...
143 
144  // calcule la bbox de la sequence de triangles
145  BBox box;
146  for(unsigned int i= begin; i < end; i++)
147  {
148  bbox_insert(box, triangles[i].a);
149  bbox_insert(box, triangles[i].b);
150  bbox_insert(box, triangles[i].c);
151  }
152 
153  // construit une feuille, s'il ne reste qu'un seul triangle
154  if(begin + 1 == end)
155  return create_leaf(box, &triangles.front() + begin);
156 
157  // trouve l'axe le plus etire de la bbox
158  // et coupe la boite en 2 sur l'axe le plus etire
159  int axis;
160  float cut;
161  Vector d= Vector(box.min, box.max);
162  if(d.x > d.y && d.x > d.z) { axis= 0; cut= (box.min.x + box.max.x) / 2; }
163  else if(d.y > d.x && d.y > d.z) { axis= 1; cut= (box.min.y + box.max.y) / 2; }
164  else { axis= 2; cut= (box.min.z + box.max.z) / 2; }
165 
166  // reparti les triangles en fonction de leur position par rapport a la coupe
167  Triangle *pmid= std::partition(&triangles.front() + begin, &triangles.front() + end, CutLess(axis, cut));
168  unsigned int mid= std::distance(&triangles.front(), pmid);
169 
170  if(mid == begin || mid == end)
171  // force un decoupage en 2 groupes, si la partition a echouee
172  mid= (begin + end) / 2;
173 
174  return create_node(box,
175  build_node(triangles, begin, mid),
176  build_node(triangles, mid, end));
177 }
178 
179 
180 BVH make_bvh( const std::vector<vec3>& positions )
181 {
182  BVH bvh;
183 
184  // construit la liste de triangles
185  bvh.triangles.reserve(positions.size() / 3);
186  for(unsigned int i= 0; i +2 < positions.size(); i= i + 3)
187  bvh.triangles.push_back( Triangle( Point(positions[i]), Point(positions[i +1]), Point(positions[i +2])) );
188 
189  printf("building bvh %d triangles...\n", (int) bvh.triangles.size());
190 
191  // construit l'arbre
192  bvh.root= build_node(bvh.triangles, 0, bvh.triangles.size());
193  bvh.box= bvh.root->box;
194 
195  return bvh;
196 }
197 
198 unsigned long int box_n= 0;
199 unsigned long int tri_n= 0;
200 
201 
202 void node_intersect( const BVHNode *node, const Ray& ray, Hit& hit )
203 {
204  if(node == nullptr) return;
205 
206  if(node->triangle)
207  {
208  //~ tri_n++;
209  // feuille, intersection avec le triangle
210  float t, u, v;
211  if(intersect(node->triangle, ray, hit.t, t, u, v))
212  {
213  // mise a jour de l'intersection courante
214  hit.t= t;
215  hit.p= ray.origin + t * ray.direction;
216 
217  #if 1
218  // calcule la normale geometrique du triangle
219  Vector ab= Vector(node->triangle->a, node->triangle->b);
220  Vector ac= Vector(node->triangle->a, node->triangle->c);
221  hit.n= normalize(cross(ab, ac));
222  #else
223  // interpole la normale au point d'intersection, necessite des normales par sommet
224  #endif
225 
226  unsigned int id= (unsigned long int) node->triangle & 255;
227  hit.color= Color(1.f - (id % 16) / 15.f, (id % 16) / 15.f, id / 255.f);
228  hit.hit= true;
229  }
230  }
231  else
232  {
233  //~ box_n++;
234  // noeud interne
235  if(intersect(node->box, ray, hit.t))
236  {
237  // visite les 2 fils
238  node_intersect(node->left, ray, hit);
239  node_intersect(node->right, ray, hit);
240  }
241  }
242 }
243 
244 bool intersect( const BVH& bvh, const Ray& ray, const float tmax, Hit& hit )
245 {
246  hit.color= Black();
247  hit.t= tmax;
248  hit.hit= false;
249  node_intersect(bvh.root, ray, hit);
250  return hit.hit;
251 }
252 
representation d'un triangle.
Definition: bvh.h:11
void end(Widgets &w)
termine la description des elements de l'interface graphique.
Definition: widgets.cpp:363
Point p
position
Definition: ray.h:37
Definition: bbox.h:10
Vector normalize(const Vector &v)
renvoie un vecteur unitaire / longueur == 1.
Definition: vec.cpp:79
Color Black()
utilitaire. renvoie une couleur noire.
Definition: color.cpp:5
representation d'une couleur (rgba) transparente ou opaque.
Definition: color.h:13
Definition: bvh.cpp:117
Color color
couleur de la surface au point d'intersection
Definition: ray.h:39
representation du bvh.
Definition: bvh.h:39
representation d'un vecteur 3d.
Definition: vec.h:42
float dot(const Vector &u, const Vector &v)
renvoie le produit scalaire de 2 vecteurs.
Definition: vec.cpp:93
float t
abscisse : p= o + t.d
Definition: ray.h:40
representation d'un point d'intersection.
Definition: ray.h:33
void begin(Widgets &w)
debut de la description des elements de l'interface graphique.
Definition: widgets.cpp:29
void printf(Text &text, const int px, const int py, const char *format,...)
affiche un texte a la position x, y. meme utilisation que printf().
Definition: text.cpp:140
Vector n
normale de la surface au point d'intersection
Definition: ray.h:38
representation d'un rayon.
Definition: ray.h:12
representation d'un noeud du bvh.
Definition: bvh.h:25
float distance(const Point &a, const Point &b)
renvoie la distance etre 2 points.
Definition: vec.cpp:7
representation d'un point 3d.
Definition: vec.h:19
Vector cross(const Vector &u, const Vector &v)
renvoie le produit vectoriel de 2 vecteurs.
Definition: vec.cpp:85