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