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   = NULL;
  54         b->primitives.end     = NULL;
  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] = NULL;
  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 (!isfinite(bb[0]) || !isfinite(bb[1]) || !isfinite(bb[2]))
 113                 return;
 114         if (!isfinite(bb[3]) || !isfinite(bb[4]) || !isfinite(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] = NULL;
 182                         tmp->sorted_end[i] = NULL;
 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;
 335                 boffset = size / 2;
 336
 337                 SweepCost *sweep = (SweepCost *)MEM_mallocN(sizeof(SweepCost) * size, "RTBuilder.HeuristicSweep");
 338
 339                 for (int axis = 0; axis < 3; axis++) {
 340                         SweepCost sweep_left;
 341
 342                         RTBuilder::Object **obj = b->sorted_begin[axis];
 343
 344 //                      float right_cost = 0;
 345                         for (int i = size - 1; i >= 0; i--) {
 346                                 if (i == size - 1) {
 347                                         copy_v3_v3(sweep[i].bb, obj[i]->bb);
 348                                         copy_v3_v3(sweep[i].bb + 3, obj[i]->bb + 3);
 349                                         sweep[i].cost = obj[i]->cost;
 350                                 }
 351                                 else {
 352                                         sweep[i].bb[0] = min_ff(obj[i]->bb[0], sweep[i + 1].bb[0]);
 353                                         sweep[i].bb[1] = min_ff(obj[i]->bb[1], sweep[i + 1].bb[1]);
 354                                         sweep[i].bb[2] = min_ff(obj[i]->bb[2], sweep[i + 1].bb[2]);
 355                                         sweep[i].bb[3] = max_ff(obj[i]->bb[3], sweep[i + 1].bb[3]);
 356                                         sweep[i].bb[4] = max_ff(obj[i]->bb[4], sweep[i + 1].bb[4]);
 357                                         sweep[i].bb[5] = max_ff(obj[i]->bb[5], sweep[i + 1].bb[5]);
 358                                         sweep[i].cost  = obj[i]->cost + sweep[i + 1].cost;
 359                                 }
 360 //                              right_cost += obj[i]->cost;
 361                         }
 362
 363                         sweep_left.bb[0] = obj[0]->bb[0];
 364                         sweep_left.bb[1] = obj[0]->bb[1];
 365                         sweep_left.bb[2] = obj[0]->bb[2];
 366                         sweep_left.bb[3] = obj[0]->bb[3];
 367                         sweep_left.bb[4] = obj[0]->bb[4];
 368                         sweep_left.bb[5] = obj[0]->bb[5];
 369                         sweep_left.cost  = obj[0]->cost;
 370
 371 //                      right_cost -= obj[0]->cost;     if (right_cost < 0) right_cost = 0;
 372
 373                         for (int i = 1; i < size; i++) {
 374                                 //Worst case heuristic (cost of each child is linear)
 375                                 float hcost, left_side, right_side;
 376
 377                                 // not using log seems to have no impact on raytracing perf, but
 378                                 // makes tree construction quicker, left out for now to test (brecht)
 379                                 // left_side  = bb_area(sweep_left.bb, sweep_left.bb + 3) * (sweep_left.cost + logf((float)i));
 380                                 // right_side = bb_area(sweep[i].bb,   sweep[i].bb   + 3) * (sweep[i].cost   + logf((float)size - i));
 381                                 left_side = bb_area(sweep_left.bb, sweep_left.bb + 3) * (sweep_left.cost);
 382                                 right_side = bb_area(sweep[i].bb, sweep[i].bb + 3) * (sweep[i].cost);
 383                                 hcost = left_side + right_side;
 384
 385                                 assert(left_side >= 0);
 386                                 assert(right_side >= 0);
 387
 388                                 if (left_side > bcost) break;   //No way we can find a better heuristic in this axis
 389
 390                                 assert(hcost >= 0);
 391                                 // this makes sure the tree built is the same whatever is the order of the sorting axis
 392                                 if (hcost < bcost || (hcost == bcost && axis < baxis)) {
 393                                         bcost = hcost;
 394                                         baxis = axis;
 395                                         boffset = i;
 396                                 }
 397                                 DO_MIN(obj[i]->bb,   sweep_left.bb);
 398                                 DO_MAX(obj[i]->bb + 3, sweep_left.bb + 3);
 399
 400                                 sweep_left.cost += obj[i]->cost;
 401 //                              right_cost -= obj[i]->cost; if (right_cost < 0) right_cost = 0;
 402                         }
 403
 404                         //assert(baxis >= 0 && baxis < 3);
 405                         if (!(baxis >= 0 && baxis < 3))
 406                                 baxis = 0;
 407                 }
 408
 409
 410                 MEM_freeN(sweep);
 411         }
 412         else if (size == 2) {
 413                 baxis = 0;
 414                 boffset = 1;
 415         }
 416         else if (size == 1) {
 417                 b->child_offset[0] = 0;
 418                 b->child_offset[1] = 1;
 419                 return 1;
 420         }
 421
 422         b->child_offset[0] = 0;
 423         b->child_offset[1] = boffset;
 424         b->child_offset[2] = size;
 425
 426
 427         /* Adjust sorted arrays for childs */
 428         for (int i = 0; i < boffset; i++) b->sorted_begin[baxis][i]->selected = true;
 429         for (int i = boffset; i < size; i++) b->sorted_begin[baxis][i]->selected = false;
 430         for (int i = 0; i < 3; i++)
 431                 std::stable_partition(b->sorted_begin[i], b->sorted_end[i], selected_node);
 432
 433         return nchilds;
 434 }
 435
 436 /*
 437  * Helper code
 438  * PARTITION code / used on mean-split
 439  * basically this a std::nth_element (like on C++ STL algorithm)
 440  */
 441 #if 0
 442 static void split_leafs(RTBuilder *b, int *nth, int partitions, int split_axis)
 443 {
 444         int i;
 445         b->split_axis = split_axis;
 446
 447         for (i = 0; i < partitions - 1; i++)
 448         {
 449                 assert(nth[i] < nth[i + 1] && nth[i + 1] < nth[partitions]);
 450
 451                 if (split_axis == 0) std::nth_element(b, nth[i],  nth[i + 1], nth[partitions], obj_bb_compare<RTBuilder::Object, 0>);
 452                 if (split_axis == 1) std::nth_element(b, nth[i],  nth[i + 1], nth[partitions], obj_bb_compare<RTBuilder::Object, 1>);
 453                 if (split_axis == 2) std::nth_element(b, nth[i],  nth[i + 1], nth[partitions], obj_bb_compare<RTBuilder::Object, 2>);
 454         }
 455 }
 456 #endif
 457
 458 /*
 459  * Bounding Box utils
 460  */
 461 float bb_volume(const float min[3], const float max[3])
 462 {
 463         return (max[0] - min[0]) * (max[1] - min[1]) * (max[2] - min[2]);
 464 }
 465
 466 float bb_area(const float min[3], const float max[3])
 467 {
 468         float sub[3], a;
 469         sub[0] = max[0] - min[0];
 470         sub[1] = max[1] - min[1];
 471         sub[2] = max[2] - min[2];
 472
 473         a = (sub[0] * sub[1] + sub[0] * sub[2] + sub[1] * sub[2]) * 2.0f;
 474         /* used to have an assert() here on negative results
 475          * however, in this case its likely some overflow or ffast math error.
 476          * so just return 0.0f instead. */
 477         return a < 0.0f ? 0.0f : a;
 478 }
 479
 480 int bb_largest_axis(const float min[3], const float max[3])
 481 {
 482         float sub[3];
 483
 484         sub[0] = max[0] - min[0];
 485         sub[1] = max[1] - min[1];
 486         sub[2] = max[2] - min[2];
 487         if (sub[0] > sub[1]) {
 488                 if (sub[0] > sub[2])
 489                         return 0;
 490                 else
 491                         return 2;
 492         }
 493         else {
 494                 if (sub[1] > sub[2])
 495                         return 1;
 496                 else
 497                         return 2;
 498         }
 499 }
 500
 501 /* only returns 0 if merging inner and outerbox would create a box larger than outer box */
 502 int bb_fits_inside(const float outer_min[3], const float outer_max[3],
 503                    const float inner_min[3], const float inner_max[3])
 504 {
 505         int i;
 506         for (i = 0; i < 3; i++)
 507                 if (outer_min[i] > inner_min[i]) return 0;
 508
 509         for (i = 0; i < 3; i++)
 510                 if (outer_max[i] < inner_max[i]) return 0;
 511
 512         return 1;
 513 }