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