51f89784674a5851ec41e7e952f6aad742bfc2af
[blender.git] / source / blender / render / intern / raytrace / rayobject_rtbuild.cpp
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2009 Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): AndrĂ© Pinto.
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/render/intern/raytrace/rayobject_rtbuild.cpp
29  *  \ingroup render
30  */
31
32
33 #include <assert.h>
34 #include <stdlib.h>
35 #include <algorithm>
36
37 #if __cplusplus >= 201103L
38 #include <cmath>
39 using std::isfinite;
40 #else
41 #include <math.h>
42 #endif
43
44 #include "rayobject_rtbuild.h"
45
46 #include "MEM_guardedalloc.h"
47
48 #include "BLI_math.h"
49 #include "BLI_utildefines.h"
50
51 static bool selected_node(RTBuilder::Object *node)
52 {
53         return node->selected;
54 }
55
56 static void rtbuild_init(RTBuilder *b)
57 {
58         b->split_axis = -1;
59         b->primitives.begin   = NULL;
60         b->primitives.end     = NULL;
61         b->primitives.maxsize = 0;
62         b->depth = 0;
63
64         for (int i = 0; i < RTBUILD_MAX_CHILDS; i++)
65                 b->child_offset[i] = 0;
66
67         for (int i = 0; i < 3; i++)
68                 b->sorted_begin[i] = b->sorted_end[i] = NULL;
69
70         INIT_MINMAX(b->bb, b->bb + 3);
71 }
72
73 RTBuilder *rtbuild_create(int size)
74 {
75         RTBuilder *builder  = (RTBuilder *) MEM_mallocN(sizeof(RTBuilder), "RTBuilder");
76         RTBuilder::Object *memblock = (RTBuilder::Object *)MEM_mallocN(sizeof(RTBuilder::Object) * size, "RTBuilder.objects");
77
78
79         rtbuild_init(builder);
80
81         builder->primitives.begin = builder->primitives.end = memblock;
82         builder->primitives.maxsize = size;
83
84         for (int i = 0; i < 3; i++) {
85                 builder->sorted_begin[i] = (RTBuilder::Object **)MEM_mallocN(sizeof(RTBuilder::Object *) * size, "RTBuilder.sorted_objects");
86                 builder->sorted_end[i]   = builder->sorted_begin[i];
87         }
88
89
90         return builder;
91 }
92
93 void rtbuild_free(RTBuilder *b)
94 {
95         if (b->primitives.begin) MEM_freeN(b->primitives.begin);
96
97         for (int i = 0; i < 3; i++)
98                 if (b->sorted_begin[i])
99                         MEM_freeN(b->sorted_begin[i]);
100
101         MEM_freeN(b);
102 }
103
104 void rtbuild_add(RTBuilder *b, RayObject *o)
105 {
106         float bb[6];
107
108         assert(b->primitives.begin + b->primitives.maxsize != b->primitives.end);
109
110         INIT_MINMAX(bb, bb + 3);
111         RE_rayobject_merge_bb(o, bb, bb + 3);
112
113         /* skip objects with invalid bounding boxes, nan causes DO_MINMAX
114          * to do nothing, so we get these invalid values. this shouldn't
115          * happen usually, but bugs earlier in the pipeline can cause it. */
116         if (bb[0] > bb[3] || bb[1] > bb[4] || bb[2] > bb[5])
117                 return;
118         /* skip objects with inf bounding boxes */
119         if (!isfinite(bb[0]) || !isfinite(bb[1]) || !isfinite(bb[2]))
120                 return;
121         if (!isfinite(bb[3]) || !isfinite(bb[4]) || !isfinite(bb[5]))
122                 return;
123         /* skip objects with zero bounding box, they are of no use, and
124          * will give problems in rtbuild_heuristic_object_split later */
125         if (bb[0] == bb[3] && bb[1] == bb[4] && bb[2] == bb[5])
126                 return;
127
128         copy_v3_v3(b->primitives.end->bb, bb);
129         copy_v3_v3(b->primitives.end->bb + 3, bb + 3);
130         b->primitives.end->obj = o;
131         b->primitives.end->cost = RE_rayobject_cost(o);
132
133         for (int i = 0; i < 3; i++) {
134                 *(b->sorted_end[i]) = b->primitives.end;
135                 b->sorted_end[i]++;
136         }
137         b->primitives.end++;
138 }
139
140 int rtbuild_size(RTBuilder *b)
141 {
142         return b->sorted_end[0] - b->sorted_begin[0];
143 }
144
145
146 template<class Obj, int Axis>
147 static bool obj_bb_compare(const Obj &a, const Obj &b)
148 {
149         if (a->bb[Axis] != b->bb[Axis])
150                 return a->bb[Axis] < b->bb[Axis];
151         return a->obj < b->obj;
152 }
153
154 template<class Item>
155 static void object_sort(Item *begin, Item *end, int axis)
156 {
157         if (axis == 0) return std::sort(begin, end, obj_bb_compare<Item, 0> );
158         if (axis == 1) return std::sort(begin, end, obj_bb_compare<Item, 1> );
159         if (axis == 2) return std::sort(begin, end, obj_bb_compare<Item, 2> );
160         assert(false);
161 }
162
163 void rtbuild_done(RTBuilder *b, RayObjectControl *ctrl)
164 {
165         for (int i = 0; i < 3; i++) {
166                 if (b->sorted_begin[i]) {
167                         if (RE_rayobjectcontrol_test_break(ctrl)) break;
168                         object_sort(b->sorted_begin[i], b->sorted_end[i], i);
169                 }
170         }
171 }
172
173 RayObject *rtbuild_get_primitive(RTBuilder *b, int index)
174 {
175         return b->sorted_begin[0][index]->obj;
176 }
177
178 RTBuilder *rtbuild_get_child(RTBuilder *b, int child, RTBuilder *tmp)
179 {
180         rtbuild_init(tmp);
181
182         tmp->depth = b->depth + 1;
183
184         for (int i = 0; i < 3; i++)
185                 if (b->sorted_begin[i]) {
186                         tmp->sorted_begin[i] = b->sorted_begin[i] +  b->child_offset[child];
187                         tmp->sorted_end[i] = b->sorted_begin[i] +  b->child_offset[child + 1];
188                 }
189                 else {
190                         tmp->sorted_begin[i] = NULL;
191                         tmp->sorted_end[i] = NULL;
192                 }
193
194         return tmp;
195 }
196
197 static void rtbuild_calc_bb(RTBuilder *b)
198 {
199         if (b->bb[0] == 1.0e30f) {
200                 for (RTBuilder::Object **index = b->sorted_begin[0]; index != b->sorted_end[0]; index++)
201                         RE_rayobject_merge_bb( (*index)->obj, b->bb, b->bb + 3);
202         }
203 }
204
205 void rtbuild_merge_bb(RTBuilder *b, float min[3], float max[3])
206 {
207         rtbuild_calc_bb(b);
208         DO_MIN(b->bb, min);
209         DO_MAX(b->bb + 3, max);
210 }
211
212 #if 0
213 int rtbuild_get_largest_axis(RTBuilder *b)
214 {
215         rtbuild_calc_bb(b);
216         return bb_largest_axis(b->bb, b->bb + 3);
217 }
218
219 //Left balanced tree
220 int rtbuild_mean_split(RTBuilder *b, int nchilds, int axis)
221 {
222         int i;
223         int mleafs_per_child, Mleafs_per_child;
224         int tot_leafs  = rtbuild_size(b);
225         int missing_leafs;
226
227         long long s;
228
229         assert(nchilds <= RTBUILD_MAX_CHILDS);
230
231         //TODO optimize calc of leafs_per_child
232         for (s = nchilds; s < tot_leafs; s *= nchilds) ;
233         Mleafs_per_child = s / nchilds;
234         mleafs_per_child = Mleafs_per_child / nchilds;
235
236         //split min leafs per child
237         b->child_offset[0] = 0;
238         for (i = 1; i <= nchilds; i++)
239                 b->child_offset[i] = mleafs_per_child;
240
241         //split remaining leafs
242         missing_leafs = tot_leafs - mleafs_per_child * nchilds;
243         for (i = 1; i <= nchilds; i++)
244         {
245                 if (missing_leafs > Mleafs_per_child - mleafs_per_child)
246                 {
247                         b->child_offset[i] += Mleafs_per_child - mleafs_per_child;
248                         missing_leafs -= Mleafs_per_child - mleafs_per_child;
249                 }
250                 else {
251                         b->child_offset[i] += missing_leafs;
252                         missing_leafs = 0;
253                         break;
254                 }
255         }
256
257         //adjust for accumulative offsets
258         for (i = 1; i <= nchilds; i++)
259                 b->child_offset[i] += b->child_offset[i - 1];
260
261         //Count created childs
262         for (i = nchilds; b->child_offset[i] == b->child_offset[i - 1]; i--) ;
263         split_leafs(b, b->child_offset, i, axis);
264
265         assert(b->child_offset[0] == 0 && b->child_offset[i] == tot_leafs);
266         return i;
267 }
268
269
270 int rtbuild_mean_split_largest_axis(RTBuilder *b, int nchilds)
271 {
272         int axis = rtbuild_get_largest_axis(b);
273         return rtbuild_mean_split(b, nchilds, axis);
274 }
275 #endif
276
277 /*
278  * "separators" is an array of dim NCHILDS-1
279  * and indicates where to cut the childs
280  */
281 #if 0
282 int rtbuild_median_split(RTBuilder *b, float *separators, int nchilds, int axis)
283 {
284         int size = rtbuild_size(b);
285
286         assert(nchilds <= RTBUILD_MAX_CHILDS);
287         if (size <= nchilds)
288         {
289                 return rtbuild_mean_split(b, nchilds, axis);
290         }
291         else {
292                 int i;
293
294                 b->split_axis = axis;
295
296                 //Calculate child offsets
297                 b->child_offset[0] = 0;
298                 for (i = 0; i < nchilds - 1; i++)
299                         b->child_offset[i + 1] = split_leafs_by_plane(b, b->child_offset[i], size, separators[i]);
300                 b->child_offset[nchilds] = size;
301
302                 for (i = 0; i < nchilds; i++)
303                         if (b->child_offset[i + 1] - b->child_offset[i] == size)
304                                 return rtbuild_mean_split(b, nchilds, axis);
305
306                 return nchilds;
307         }
308 }
309
310 int rtbuild_median_split_largest_axis(RTBuilder *b, int nchilds)
311 {
312         int la, i;
313         float separators[RTBUILD_MAX_CHILDS];
314
315         rtbuild_calc_bb(b);
316
317         la = bb_largest_axis(b->bb, b->bb + 3);
318         for (i = 1; i < nchilds; i++)
319                 separators[i - 1] = (b->bb[la + 3] - b->bb[la]) * i / nchilds;
320
321         return rtbuild_median_split(b, separators, nchilds, la);
322 }
323 #endif
324
325 //Heuristics Object Splitter
326
327
328 struct SweepCost {
329         float bb[6];
330         float cost;
331 };
332
333 /* Object Surface Area Heuristic splitter */
334 int rtbuild_heuristic_object_split(RTBuilder *b, int nchilds)
335 {
336         int size = rtbuild_size(b);
337         assert(nchilds == 2);
338         assert(size > 1);
339         int baxis = -1, boffset = 0;
340
341         if (size > nchilds) {
342                 if (b->depth > RTBUILD_MAX_SAH_DEPTH) {
343                         // for degenerate cases we avoid running out of stack space
344                         // by simply splitting the children in the middle
345                         b->child_offset[0] = 0;
346                         b->child_offset[1] = (size+1)/2;
347                         b->child_offset[2] = size;
348                         return 2;
349                 }
350
351                 float bcost = FLT_MAX;
352                 baxis = -1;
353                 boffset = size / 2;
354
355                 SweepCost *sweep = (SweepCost *)MEM_mallocN(sizeof(SweepCost) * size, "RTBuilder.HeuristicSweep");
356
357                 for (int axis = 0; axis < 3; axis++) {
358                         SweepCost sweep_left;
359
360                         RTBuilder::Object **obj = b->sorted_begin[axis];
361
362 //                      float right_cost = 0;
363                         for (int i = size - 1; i >= 0; i--) {
364                                 if (i == size - 1) {
365                                         copy_v3_v3(sweep[i].bb, obj[i]->bb);
366                                         copy_v3_v3(sweep[i].bb + 3, obj[i]->bb + 3);
367                                         sweep[i].cost = obj[i]->cost;
368                                 }
369                                 else {
370                                         sweep[i].bb[0] = min_ff(obj[i]->bb[0], sweep[i + 1].bb[0]);
371                                         sweep[i].bb[1] = min_ff(obj[i]->bb[1], sweep[i + 1].bb[1]);
372                                         sweep[i].bb[2] = min_ff(obj[i]->bb[2], sweep[i + 1].bb[2]);
373                                         sweep[i].bb[3] = max_ff(obj[i]->bb[3], sweep[i + 1].bb[3]);
374                                         sweep[i].bb[4] = max_ff(obj[i]->bb[4], sweep[i + 1].bb[4]);
375                                         sweep[i].bb[5] = max_ff(obj[i]->bb[5], sweep[i + 1].bb[5]);
376                                         sweep[i].cost  = obj[i]->cost + sweep[i + 1].cost;
377                                 }
378 //                              right_cost += obj[i]->cost;
379                         }
380
381                         sweep_left.bb[0] = obj[0]->bb[0];
382                         sweep_left.bb[1] = obj[0]->bb[1];
383                         sweep_left.bb[2] = obj[0]->bb[2];
384                         sweep_left.bb[3] = obj[0]->bb[3];
385                         sweep_left.bb[4] = obj[0]->bb[4];
386                         sweep_left.bb[5] = obj[0]->bb[5];
387                         sweep_left.cost  = obj[0]->cost;
388
389 //                      right_cost -= obj[0]->cost;     if (right_cost < 0) right_cost = 0;
390
391                         for (int i = 1; i < size; i++) {
392                                 //Worst case heuristic (cost of each child is linear)
393                                 float hcost, left_side, right_side;
394
395                                 // not using log seems to have no impact on raytracing perf, but
396                                 // makes tree construction quicker, left out for now to test (brecht)
397                                 // left_side  = bb_area(sweep_left.bb, sweep_left.bb + 3) * (sweep_left.cost + logf((float)i));
398                                 // right_side = bb_area(sweep[i].bb,   sweep[i].bb   + 3) * (sweep[i].cost   + logf((float)size - i));
399                                 left_side = bb_area(sweep_left.bb, sweep_left.bb + 3) * (sweep_left.cost);
400                                 right_side = bb_area(sweep[i].bb, sweep[i].bb + 3) * (sweep[i].cost);
401                                 hcost = left_side + right_side;
402
403                                 assert(left_side >= 0);
404                                 assert(right_side >= 0);
405
406                                 if (left_side > bcost) break;   //No way we can find a better heuristic in this axis
407
408                                 assert(hcost >= 0);
409                                 // this makes sure the tree built is the same whatever is the order of the sorting axis
410                                 if (hcost < bcost || (hcost == bcost && axis < baxis)) {
411                                         bcost = hcost;
412                                         baxis = axis;
413                                         boffset = i;
414                                 }
415                                 DO_MIN(obj[i]->bb,   sweep_left.bb);
416                                 DO_MAX(obj[i]->bb + 3, sweep_left.bb + 3);
417
418                                 sweep_left.cost += obj[i]->cost;
419 //                              right_cost -= obj[i]->cost; if (right_cost < 0) right_cost = 0;
420                         }
421
422                         //assert(baxis >= 0 && baxis < 3);
423                         if (!(baxis >= 0 && baxis < 3))
424                                 baxis = 0;
425                 }
426
427
428                 MEM_freeN(sweep);
429         }
430         else if (size == 2) {
431                 baxis = 0;
432                 boffset = 1;
433         }
434         else if (size == 1) {
435                 b->child_offset[0] = 0;
436                 b->child_offset[1] = 1;
437                 return 1;
438         }
439
440         b->child_offset[0] = 0;
441         b->child_offset[1] = boffset;
442         b->child_offset[2] = size;
443
444
445         /* Adjust sorted arrays for childs */
446         for (int i = 0; i < boffset; i++) b->sorted_begin[baxis][i]->selected = true;
447         for (int i = boffset; i < size; i++) b->sorted_begin[baxis][i]->selected = false;
448         for (int i = 0; i < 3; i++)
449                 std::stable_partition(b->sorted_begin[i], b->sorted_end[i], selected_node);
450
451         return nchilds;
452 }
453
454 /*
455  * Helper code
456  * PARTITION code / used on mean-split
457  * basically this a std::nth_element (like on C++ STL algorithm)
458  */
459 #if 0
460 static void split_leafs(RTBuilder *b, int *nth, int partitions, int split_axis)
461 {
462         int i;
463         b->split_axis = split_axis;
464
465         for (i = 0; i < partitions - 1; i++)
466         {
467                 assert(nth[i] < nth[i + 1] && nth[i + 1] < nth[partitions]);
468
469                 if (split_axis == 0) std::nth_element(b, nth[i],  nth[i + 1], nth[partitions], obj_bb_compare<RTBuilder::Object, 0>);
470                 if (split_axis == 1) std::nth_element(b, nth[i],  nth[i + 1], nth[partitions], obj_bb_compare<RTBuilder::Object, 1>);
471                 if (split_axis == 2) std::nth_element(b, nth[i],  nth[i + 1], nth[partitions], obj_bb_compare<RTBuilder::Object, 2>);
472         }
473 }
474 #endif
475
476 /*
477  * Bounding Box utils
478  */
479 float bb_volume(const float min[3], const float max[3])
480 {
481         return (max[0] - min[0]) * (max[1] - min[1]) * (max[2] - min[2]);
482 }
483
484 float bb_area(const float min[3], const float max[3])
485 {
486         float sub[3], a;
487         sub[0] = max[0] - min[0];
488         sub[1] = max[1] - min[1];
489         sub[2] = max[2] - min[2];
490
491         a = (sub[0] * sub[1] + sub[0] * sub[2] + sub[1] * sub[2]) * 2.0f;
492         /* used to have an assert() here on negative results
493          * however, in this case its likely some overflow or ffast math error.
494          * so just return 0.0f instead. */
495         return a < 0.0f ? 0.0f : a;
496 }
497
498 int bb_largest_axis(const float min[3], const float max[3])
499 {
500         float sub[3];
501
502         sub[0] = max[0] - min[0];
503         sub[1] = max[1] - min[1];
504         sub[2] = max[2] - min[2];
505         if (sub[0] > sub[1]) {
506                 if (sub[0] > sub[2])
507                         return 0;
508                 else
509                         return 2;
510         }
511         else {
512                 if (sub[1] > sub[2])
513                         return 1;
514                 else
515                         return 2;
516         }
517 }
518
519 /* only returns 0 if merging inner and outerbox would create a box larger than outer box */
520 int bb_fits_inside(const float outer_min[3], const float outer_max[3],
521                    const float inner_min[3], const float inner_max[3])
522 {
523         int i;
524         for (i = 0; i < 3; i++)
525                 if (outer_min[i] > inner_min[i]) return 0;
526
527         for (i = 0; i < 3; i++)
528                 if (outer_max[i] < inner_max[i]) return 0;
529
530         return 1;
531 }