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