update for MingW/CMake
[blender-staging.git] / source / blender / render / intern / raytrace / reorganize.h
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 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 <stdio.h>
30 #include <algorithm>
31 #include <math.h>
32 #include <vector>
33 #include <queue>
34
35 #include "BKE_global.h"
36
37 #ifdef _WIN32
38 #undef INFINITY
39 #define INFINITY FLT_MAX // in mingw math.h: (1.0F/0.0F). This generates compile error, though.
40 #endif
41
42 extern int tot_pushup;
43 extern int tot_pushdown;
44
45 #if !defined(INFINITY) && defined(HUGE_VAL)
46 #define INFINITY HUGE_VAL
47 #endif
48
49 template<class Node>
50 bool node_fits_inside(Node *a, Node *b)
51 {
52         return bb_fits_inside(b->bb, b->bb+3, a->bb, a->bb+3);
53 }
54
55 template<class Node>
56 void reorganize_find_fittest_parent(Node *tree, Node *node, std::pair<float,Node*> &cost)
57 {
58         std::queue<Node*> q;
59         q.push(tree);
60         
61         while(!q.empty())
62         {
63                 Node *parent = q.front();
64                 q.pop();
65                 
66                 if(parent == node) continue;
67                 if(node_fits_inside(node, parent) && RE_rayobject_isAligned(parent->child) )
68                 {
69                         float pcost = bb_area(parent->bb, parent->bb+3);
70                         cost = std::min( cost, std::make_pair(pcost,parent) );
71                         for(Node *child = parent->child; child; child = child->sibling)
72                                 q.push(child);                  
73                 }
74         }
75 }
76
77 static int tot_moves = 0;
78 template<class Node>
79 void reorganize(Node *root)
80 {
81         std::queue<Node*> q;
82
83         q.push(root);
84         while(!q.empty())
85         {
86                 Node * node = q.front();
87                 q.pop();
88                 
89                 if( RE_rayobject_isAligned(node->child) )
90                 {
91                         for(Node **prev = &node->child; *prev; )
92                         {
93                                 assert( RE_rayobject_isAligned(*prev) ); 
94                                 q.push(*prev);
95
96                                 std::pair<float,Node*> best(FLT_MAX, root);
97                                 reorganize_find_fittest_parent( root, *prev, best );
98
99                                 if(best.second == node)
100                                 {
101                                         //Already inside the fitnest BB
102                                         prev = &(*prev)->sibling;
103                                 }
104                                 else
105                                 {
106                                         Node *tmp = *prev;
107                                         *prev = (*prev)->sibling;
108                                         
109                                         tmp->sibling =  best.second->child;
110                                         best.second->child = tmp;
111                                         
112                                         tot_moves++;
113                                 }
114                         
115                         
116                         }
117                 }
118                 if(node != root)
119                 {
120                 }
121         }
122 }
123
124 /*
125  * Prunes useless nodes from trees:
126  *  erases nodes with total ammount of primitives = 0
127  *  prunes nodes with only one child (except if that child is a primitive)
128  */
129 template<class Node>
130 void remove_useless(Node *node, Node **new_node)
131 {
132         if( RE_rayobject_isAligned(node->child) )
133         {
134
135                 for(Node **prev = &node->child; *prev; )
136                 {
137                         Node *next = (*prev)->sibling;
138                         remove_useless(*prev, prev);
139                         if(*prev == 0)
140                                 *prev = next;
141                         else
142                         {
143                                 (*prev)->sibling = next;
144                                 prev = &((*prev)->sibling);
145                         }
146                 }                       
147         }
148         if(node->child)
149         {
150                 if(RE_rayobject_isAligned(node->child) && node->child->sibling == 0)
151                         *new_node = node->child;
152         }
153         else if(node->child == 0)
154                 *new_node = 0;  
155 }
156
157 /*
158  * Minimizes expected number of BBtest by colapsing nodes
159  * it uses surface area heuristic for determining whether a node should be colapsed
160  */
161 template<class Node>
162 void pushup(Node *parent)
163 {
164         if(is_leaf(parent)) return;
165         
166         float p_area = bb_area(parent->bb, parent->bb+3);
167         Node **prev = &parent->child;
168         for(Node *child = parent->child; RE_rayobject_isAligned(child) && child; )
169         {
170                 float c_area = bb_area(child->bb, child->bb+3) ;
171                 int nchilds = count_childs(child);
172                 float original_cost = ((p_area != 0.0f)? (c_area / p_area)*nchilds: 1.0f) + 1;
173                 float flatten_cost = nchilds;
174                 if(flatten_cost < original_cost && nchilds >= 2)
175                 {
176                         append_sibling(child, child->child);
177                         child = child->sibling;
178                         *prev = child;
179
180 //                      *prev = child->child;
181 //                      append_sibling( *prev, child->sibling );
182 //                      child = *prev;
183                         tot_pushup++;
184                 }
185                 else
186                 {
187                         *prev = child;
188                         prev = &(*prev)->sibling;
189                         child = *prev;
190                 }               
191         }
192         
193         for(Node *child = parent->child; RE_rayobject_isAligned(child) && child; child = child->sibling)
194                 pushup(child);
195 }
196
197 /*
198  * try to optimize number of childs to be a multiple of SSize
199  */
200 template<class Node, int SSize>
201 void pushup_simd(Node *parent)
202 {
203         if(is_leaf(parent)) return;
204         
205         int n = count_childs(parent);
206                 
207         Node **prev = &parent->child;
208         for(Node *child = parent->child; RE_rayobject_isAligned(child) && child; )
209         {
210                 int cn = count_childs(child);
211                 if(cn-1 <= (SSize - (n%SSize) ) % SSize && RE_rayobject_isAligned(child->child) )
212                 {
213                         n += (cn - 1);
214                         append_sibling(child, child->child);
215                         child = child->sibling;
216                         *prev = child;  
217                 }
218                 else
219                 {
220                         *prev = child;
221                         prev = &(*prev)->sibling;
222                         child = *prev;
223                 }               
224         }
225                 
226         for(Node *child = parent->child; RE_rayobject_isAligned(child) && child; child = child->sibling)
227                 pushup_simd<Node,SSize>(child);
228 }
229
230
231 /*
232  * Pushdown
233  *      makes sure no child fits inside any of its sibling
234  */
235 template<class Node>
236 void pushdown(Node *parent)
237 {
238         Node **s_child = &parent->child;
239         Node * child = parent->child;
240         
241         while(child && RE_rayobject_isAligned(child))
242         {
243                 Node *next = child->sibling;
244                 Node **next_s_child = &child->sibling;
245                 
246                 //assert(bb_fits_inside(parent->bb, parent->bb+3, child->bb, child->bb+3));
247                 
248                 for(Node *i = parent->child; RE_rayobject_isAligned(i) && i; i = i->sibling)
249                 if(child != i && bb_fits_inside(i->bb, i->bb+3, child->bb, child->bb+3) && RE_rayobject_isAligned(i->child))
250                 {
251 //                      todo optimize (should the one with the smallest area?)
252 //                      float ia = bb_area(i->bb, i->bb+3)
253 //                      if(child->i)
254                         *s_child = child->sibling;
255                         child->sibling = i->child;
256                         i->child = child;
257                         next_s_child = s_child;
258                         
259                         tot_pushdown++;
260                         break;
261                 }
262                 child = next;
263                 s_child = next_s_child;
264         }
265         
266         for(Node *i = parent->child; RE_rayobject_isAligned(i) && i; i = i->sibling)
267                 pushdown( i );  
268 }
269
270
271 /*
272  * BVH refit
273  * readjust nodes BB (useful if nodes childs where modified)
274  */
275 template<class Node>
276 float bvh_refit(Node *node)
277 {
278         if(is_leaf(node)) return 0;     
279         if(is_leaf(node->child)) return 0;
280         
281         float total = 0;
282         
283         for(Node *child = node->child; child; child = child->sibling)
284                 total += bvh_refit(child);
285                 
286         float old_area = bb_area(node->bb, node->bb+3);
287         INIT_MINMAX(node->bb, node->bb+3);
288         for(Node *child = node->child; child; child = child->sibling)
289         {
290                 DO_MIN(child->bb, node->bb);
291                 DO_MAX(child->bb+3, node->bb+3);
292         }
293         total += old_area - bb_area(node->bb, node->bb+3);
294         return total;
295 }
296
297
298 /*
299  * this finds the best way to packing a tree according to a given test cost function
300  * with the purpose to reduce the expected cost (eg.: number of BB tests).
301  */
302 #include <vector>
303 #define MAX_CUT_SIZE    16
304 #define MAX_OPTIMIZE_CHILDS     MAX_CUT_SIZE
305
306 struct OVBVHNode
307 {
308         float   bb[6];
309
310         OVBVHNode *child;
311         OVBVHNode *sibling;
312         
313         /*
314          * Returns min cost to represent the subtree starting at the given node,
315          * allowing it to have a given cutsize
316          */
317         float cut_cost[MAX_CUT_SIZE];
318         float get_cost(int cutsize)
319         {
320                 return cut_cost[cutsize-1];
321         }
322         
323         /*
324          * This saves the cut size of this child, when parent is reaching
325          * its minimum cut with the given cut size
326          */
327         int cut_size[MAX_CUT_SIZE];
328         int get_cut_size(int parent_cut_size)
329         {
330                 return cut_size[parent_cut_size-1];
331         }
332         
333         /*
334          * Reorganize the node based on calculated cut costs
335          */      
336         int best_cutsize;
337         void set_cut(int cutsize, OVBVHNode ***cut)
338         {
339                 if(cutsize == 1)
340                 {
341                         **cut = this;
342                          *cut = &(**cut)->sibling;
343                 }
344                 else
345                 {
346                         if(cutsize > MAX_CUT_SIZE)
347                         {
348                                 for(OVBVHNode *child = this->child; child && RE_rayobject_isAligned(child); child = child->sibling)
349                                 {
350                                         child->set_cut( 1, cut );
351                                         cutsize--;
352                                 }
353                                 assert(cutsize == 0);
354                         }
355                         else
356                                 for(OVBVHNode *child = this->child; child && RE_rayobject_isAligned(child); child = child->sibling)
357                                         child->set_cut( child->get_cut_size( cutsize ), cut );
358                 }
359         }
360
361         void optimize()
362         {
363                 if(RE_rayobject_isAligned(this->child))
364                 {
365                         //Calc new childs
366                         {
367                                 OVBVHNode **cut = &(this->child);
368                                 set_cut( best_cutsize, &cut );
369                                 *cut = NULL;
370                         }
371
372                         //Optimize new childs
373                         for(OVBVHNode *child = this->child; child && RE_rayobject_isAligned(child); child = child->sibling)
374                                 child->optimize();
375                 }               
376         }
377 };
378
379 /*
380  * Calculates an optimal SIMD packing
381  *
382  */
383 template<class Node, class TestCost>
384 struct VBVH_optimalPackSIMD
385 {
386         TestCost testcost;
387         
388         VBVH_optimalPackSIMD(TestCost testcost)
389         {
390                 this->testcost = testcost;
391         }
392         
393         /*
394          * calc best cut on a node
395          */
396         struct calc_best
397         {
398                 Node *child[MAX_OPTIMIZE_CHILDS];
399                 float child_hit_prob[MAX_OPTIMIZE_CHILDS];
400                 
401                 calc_best(Node *node)
402                 {
403                         int nchilds = 0;
404                         //Fetch childs and needed data
405                         {
406                                 float parent_area = bb_area(node->bb, node->bb+3);
407                                 for(Node *child = node->child; child && RE_rayobject_isAligned(child); child = child->sibling)
408                                 {
409                                         this->child[nchilds] = child;
410                                         this->child_hit_prob[nchilds] = (parent_area != 0.0f)? bb_area(child->bb, child->bb+3) / parent_area: 1.0f;
411                                         nchilds++;
412                                 }
413
414                                 assert(nchilds >= 2 && nchilds <= MAX_OPTIMIZE_CHILDS);
415                         }
416                         
417                         
418                         //Build DP table to find minimum cost to represent this node with a given cutsize
419                         int   bt  [MAX_OPTIMIZE_CHILDS+1][MAX_CUT_SIZE+1]; //backtrace table
420                         float cost[MAX_OPTIMIZE_CHILDS+1][MAX_CUT_SIZE+1]; //cost table (can be reduced to float[2][MAX_CUT_COST])
421                         
422                         for(int i=0; i<=nchilds; i++)
423                         for(int j=0; j<=MAX_CUT_SIZE; j++)
424                                 cost[i][j] = INFINITY;
425
426                         cost[0][0] = 0;
427                         
428                         for(int i=1; i<=nchilds; i++)
429                         for(int size=i-1; size/*+(nchilds-i)*/<=MAX_CUT_SIZE; size++)
430                         for(int cut=1; cut+size/*+(nchilds-i)*/<=MAX_CUT_SIZE; cut++)
431                         {
432                                 float new_cost = cost[i-1][size] + child_hit_prob[i-1]*child[i-1]->get_cost(cut);
433                                 if(new_cost < cost[i][size+cut])
434                                 {
435                                         cost[i][size+cut] = new_cost;
436                                         bt[i][size+cut] = cut;
437                                 }
438                         }
439                         
440                         //Save the ways to archieve the minimum cost with a given cutsize
441                         for(int i = nchilds; i <= MAX_CUT_SIZE; i++)
442                         {
443                                 node->cut_cost[i-1] = cost[nchilds][i];
444                                 if(cost[nchilds][i] < INFINITY)
445                                 {
446                                         int current_size = i;
447                                         for(int j=nchilds; j>0; j--)
448                                         {
449                                                 child[j-1]->cut_size[i-1] = bt[j][current_size];
450                                                 current_size -= bt[j][current_size];
451                                         }
452                                 }
453                         }                       
454                 }
455         };
456         
457         void calc_costs(Node *node)
458         {
459                 
460                 if( RE_rayobject_isAligned(node->child) )
461                 {
462                         int nchilds = 0;
463                         for(Node *child = node->child; child && RE_rayobject_isAligned(child); child = child->sibling)
464                         {
465                                 calc_costs(child);
466                                 nchilds++;
467                         }
468
469                         for(int i=0; i<MAX_CUT_SIZE; i++)
470                                 node->cut_cost[i] = INFINITY;
471
472                         //We are not allowed to look on nodes with with so many childs
473                         if(nchilds > MAX_CUT_SIZE)
474                         {
475                                 float cost = 0;
476
477                                 float parent_area = bb_area(node->bb, node->bb+3);
478                                 for(Node *child = node->child; child && RE_rayobject_isAligned(child); child = child->sibling)
479                                 {
480                                         cost += ((parent_area != 0.0f)? ( bb_area(child->bb, child->bb+3) / parent_area ): 1.0f) * child->get_cost(1);
481                                 }
482                                 
483                                 cost += testcost(nchilds);
484                                 node->cut_cost[0] = cost;
485                                 node->best_cutsize = nchilds;
486                         }
487                         else
488                         {
489                                 calc_best calc(node);
490                 
491                                 //calc expected cost if we optimaly pack this node
492                                 for(int cutsize=nchilds; cutsize<=MAX_CUT_SIZE; cutsize++)
493                                 {
494                                         float m = node->get_cost(cutsize) + testcost(cutsize);
495                                         if(m < node->cut_cost[0])
496                                         {
497                                                 node->cut_cost[0] = m;
498                                                 node->best_cutsize = cutsize;
499                                         }
500                                 }
501                         }
502                         assert(node->cut_cost[0] != INFINITY);
503                 }
504                 else
505                 {
506                         node->cut_cost[0] = 1.0f;
507                         for(int i=1; i<MAX_CUT_SIZE; i++)
508                                 node->cut_cost[i] = INFINITY;
509                 }
510         }
511
512         Node *transform(Node *node)
513         {
514                 if(RE_rayobject_isAligned(node->child))
515                 {
516                         static int num = 0;
517                         bool first = false;
518                         if(num == 0) { num++; first = true; }
519                         
520                         calc_costs(node);
521                         if((G.f & G_DEBUG) && first) printf("expected cost = %f (%d)\n", node->cut_cost[0], node->best_cutsize );
522                         node->optimize();
523                 }
524                 return node;            
525         }       
526 };