diff -r -u -N -x linux2-config.py -x CVS blender-cvs/intern/guardedalloc/intern/mallocn.c blender/intern/guardedalloc/intern/mallocn.c --- blender-cvs/intern/guardedalloc/intern/mallocn.c 2006-06-16 15:27:03.000000000 -0400 +++ blender/intern/guardedalloc/intern/mallocn.c 2006-06-16 13:34:32.000000000 -0400 @@ -39,6 +39,7 @@ #include #include /* memcpy */ #include +#include /* mmap exception */ #if defined(AMIGA) || defined(__BeOS) || defined(WIN32) @@ -234,6 +235,69 @@ return NULL; } +void *MEM_reallocN(void *ptr, unsigned int len) +{ + MemHead *memh; + MemTail *memt; + MemTail memt_tmp; + static const char *name; + + /* emulate real realloc */ + if (ptr == NULL) + return (void*) MEM_mallocN(len, "realloc NULL ptr"); + + /* Recover structure pointers */ + memh = ((MemHead*) ptr) - 1; + memt = (MemTail *)(((char *) memh) + sizeof(MemHead) + memh->len); + + /* Sanity checks */ + if(sizeof(long)==8) { + if (((long) memh) & 0x7) { + MemorY_ErroR("realloc","attempt to realloc illegal pointer"); + return NULL; + } + } else { + if (((long) memh) & 0x3) { + MemorY_ErroR("realloc","attempt to realloc illegal pointer"); + return NULL; + } + } + + if ((memh->tag1 == MEMTAG1) && (memh->tag2 == MEMTAG2) && ((memh->len & 0x3) == 0)) { + if (memt->tag3 != MEMTAG3){ + MemorY_ErroR(memh->name, "end corrupt"); + name = check_memlist(memh); + if (name != 0 && name != memh->name) + MemorY_ErroR(name,"is also corrupt"); + } + } else { + name = check_memlist(memh); + if (name == 0) MemorY_ErroR("realloc","pointer not in memlist"); + else MemorY_ErroR(name,"realloc error in header"); + } + + /* Save old memt because it could get dropped in the realloc */ + memt_tmp.tag3 = memt->tag3; + memt_tmp.pad = memt->pad; + + len = (len + 3 ) & ~3; /* allocate in units of 4 */ + + memh = (MemHead *)realloc(memh, len+sizeof(MemHead)+sizeof(MemTail)); + memt = (MemTail *)(((char *) memh) + sizeof(MemHead) + len); + + /* Perform surgery on existing memt / memh */ + if(memh) { + memt = (MemTail *)(((char *) memh) + sizeof(MemHead) + len); + memt->tag3 = memt_tmp.tag3; + memt->pad = memt_tmp.pad; + mem_in_use += len - memh->len; + memh->len = len; + return ptr; + } + print_error("Realloc returns nill: len=%d in %s, total %u\n",len, memh->name, mem_in_use); + return NULL; +} + void *MEM_callocN(unsigned int len, const char *str) { MemHead *memh; diff -r -u -N -x linux2-config.py -x CVS blender-cvs/intern/guardedalloc/MEM_guardedalloc.h blender/intern/guardedalloc/MEM_guardedalloc.h --- blender-cvs/intern/guardedalloc/MEM_guardedalloc.h 2006-06-16 15:27:03.000000000 -0400 +++ blender/intern/guardedalloc/MEM_guardedalloc.h 2006-06-16 02:48:24.000000000 -0400 @@ -99,6 +99,11 @@ * */ void *MEM_mapallocN(unsigned int len, const char * str); + /** Resize previously allocated memory at address ptr to size len. + * Can be free'd with MEM_freeN as usual. + * */ + void *MEM_reallocN(void *ptr, unsigned int len); + /** Print a list of the names and sizes of all allocated memory * blocks. */ void MEM_printmemlist(void); diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/include/define.h blender/source/blender/render/intern/include/define.h --- blender-cvs/source/blender/render/intern/include/define.h 1969-12-31 19:00:00.000000000 -0500 +++ blender/source/blender/render/intern/include/define.h 2006-06-14 03:39:06.000000000 -0400 @@ -0,0 +1,247 @@ +/* D E F I N E . H + * + * @file define.h + * + * BRL-CAD + * + * Copyright (c) 2002-2006 United States Government as represented by + * the U.S. Army Research Laboratory. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2 of + * the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this file; see the file named COPYING for more + * information. + * + * D E F I N E. H + * + * Comments - + * Triangle Intersection Engine Defines + * + * Author - + * Justin L. Shumaker + * + * Source - + * The U. S. Army Research Laboratory + * Aberdeen Proving Ground, Maryland 21005-5068 USA + * + * $Id: define.h,v 1.23 2006/01/18 06:46:11 brlcad Exp $ + */ +/** @addtogroup libtie */ /** @{ */ + +#ifndef _TIE_DEFINE_H +#define _TIE_DEFINE_H + +#include + +#define TIE_SINGLE_PREC 1 /* Use Single Precision Math */ +#define TIE_TAB1 "\1\0\0\2\2\1" /* Triangle Index Table */ +#define TIE_KDTREE_NODE_MAX 4 /* Maximum number of triangles that can reside in a given node until it should be split */ +#define TIE_KDTREE_DEPTH_K1 1.6 /* K1 Depth Constant Coefficient */ +#define TIE_KDTREE_DEPTH_K2 1 /* K2 Contant */ + +#define MAX_SLICES 100 +#define MIN_SLICES 35 +#define MIN_DENSITY 0.01 +#define MIN_SPAN 0.15 +#define SCALE_COEF 1.80 + + +/* Type to use for floating precision */ +#if TIE_SINGLE_PREC +#define tfloat float +#else +#define tfloat double +#endif + + +/** + * + */ +#define math_min3(_a, _b, _c, _d) { \ + _a = _b < _c ? _b < _d ? _b : _d : _c < _d ? _c : _d; } + +/** + * + */ +#define math_max3(_a, _b, _c, _d) { \ + _a = _b > _c ? _b > _d ? _b : _d : _c > _d ? _c : _d; } + +/** + * + */ +#define math_vec_set(_a, _b, _c, _d) { \ + _a.v[0] = _b; \ + _a.v[1] = _c; \ + _a.v[2] = _d; } + +/** + * + */ +#define math_vec_min(_a, _b) { \ + _a.v[0] = _a.v[0] < _b.v[0] ? _a.v[0] : _b.v[0]; \ + _a.v[1] = _a.v[1] < _b.v[1] ? _a.v[1] : _b.v[1]; \ + _a.v[2] = _a.v[2] < _b.v[2] ? _a.v[2] : _b.v[2]; } + +/** + * + */ +#define math_vec_max(_a, _b) { \ + _a.v[0] = _a.v[0] > _b.v[0] ? _a.v[0] : _b.v[0]; \ + _a.v[1] = _a.v[1] > _b.v[1] ? _a.v[1] : _b.v[1]; \ + _a.v[2] = _a.v[2] > _b.v[2] ? _a.v[2] : _b.v[2]; } + +/** + * + */ +#define math_vec_cross(_a, _b, _c) {\ + _a.v[0] = _b.v[1]*_c.v[2] - _b.v[2]*_c.v[1]; \ + _a.v[1] = _b.v[2]*_c.v[0] - _b.v[0]*_c.v[2]; \ + _a.v[2] = _b.v[0]*_c.v[1] - _b.v[1]*_c.v[0]; } + +/** + * + */ +#define math_vec_unitize(_a) { \ + tfloat _b = 1/sqrt(_a.v[0]*_a.v[0] + _a.v[1]*_a.v[1] + _a.v[2]*_a.v[2]); \ + _a.v[0] *= _b; _a.v[1] *= _b; _a.v[2] *= _b; } + +/** + * + */ +#define math_vec_sq(_a, _b) { \ + _a.v[0] = _b.v[0] * _b.v[0]; \ + _a.v[1] = _b.v[1] * _b.v[1]; \ + _a.v[2] = _b.v[2] * _b.v[2]; } + +/** + * + */ +#define math_vec_dot(_a, _b, _c) { \ + _a = _b.v[0]*_c.v[0] + _b.v[1]*_c.v[1] + _b.v[2]*_c.v[2]; } + +/** + * + */ +#define math_vec_add(_a, _b, _c) { \ + _a.v[0] = _b.v[0] + _c.v[0]; \ + _a.v[1] = _b.v[1] + _c.v[1]; \ + _a.v[2] = _b.v[2] + _c.v[2]; } + +/** + * + */ +#define math_vec_sub(_a, _b, _c) { \ + _a.v[0] = _b.v[0] - _c.v[0]; \ + _a.v[1] = _b.v[1] - _c.v[1]; \ + _a.v[2] = _b.v[2] - _c.v[2]; } + +/** + * + */ +#define math_vec_mul(_a, _b, _c) { \ + _a.v[0] = _b.v[0] * _c.v[0]; \ + _a.v[1] = _b.v[1] * _c.v[1]; \ + _a.v[2] = _b.v[2] * _c.v[2]; } + +/** + * + */ +#define math_vec_mul_scalar(_a, _b, _c) { \ + _a.v[0] = _b.v[0] * _c; \ + _a.v[1] = _b.v[1] * _c; \ + _a.v[2] = _b.v[2] * _c; } + +/** + * + */ +#define math_vec_div(_a, _b, _c) { \ + _a.v[0] = _b.v[0] / _c.v[0]; \ + _a.v[1] = _b.v[1] / _c.v[1]; \ + _a.v[2] = _b.v[2] / _c.v[2]; } + +/** + * + */ +#define math_bbox(_a, _b, _c, _d, _e) { \ + math_min3(_a.v[0], _c.v[0], _d.v[0], _e.v[0]); \ + math_min3(_a.v[1], _c.v[1], _d.v[1], _e.v[1]); \ + math_min3(_a.v[2], _c.v[2], _d.v[2], _e.v[2]); \ + math_max3(_b.v[0], _c.v[0], _d.v[0], _e.v[0]); \ + math_max3(_b.v[1], _c.v[1], _d.v[1], _e.v[1]); \ + math_max3(_b.v[2], _c.v[2], _d.v[2], _e.v[2]); } + +/* ======================== X-tests ======================== */ +/** + * + */ +#define AXISTEST_X01(a, b, fa, fb) \ + p.v[0] = a*v0.v[1] - b*v0.v[2]; \ + p.v[2] = a*v2.v[1] - b*v2.v[2]; \ + if (p.v[0] < p.v[2]) { min = p.v[0]; max = p.v[2]; } else { min = p.v[2]; max = p.v[0]; } \ + rad = fa * half_size -> v[1] + fb * half_size -> v[2]; \ + if (min > rad || max < -rad) return 0; \ + +/** + * + */ +#define AXISTEST_X2(a, b, fa, fb) \ + p.v[0] = a*v0.v[1] - b*v0.v[2]; \ + p.v[1] = a*v1.v[1] - b*v1.v[2]; \ + if (p.v[0] < p.v[1]) { min = p.v[0]; max = p.v[1]; } else { min = p.v[1]; max = p.v[0]; } \ + rad = fa * half_size -> v[1] + fb * half_size -> v[2]; \ + if (min > rad || max < -rad) return 0; + +/* ======================== Y-tests ======================== */ +/** + * + */ +#define AXISTEST_Y02(a, b, fa, fb) \ + p.v[0] = -a*v0.v[0] + b*v0.v[2]; \ + p.v[2] = -a*v2.v[0] + b*v2.v[2]; \ + if (p.v[0] < p.v[2]) { min = p.v[0]; max = p.v[2]; } else { min = p.v[2]; max = p.v[0]; } \ + rad = fa * half_size -> v[0] + fb * half_size -> v[2]; \ + if (min > rad || max < -rad) return 0; + +/** + * + */ +#define AXISTEST_Y1(a, b, fa, fb) \ + p.v[0] = -a*v0.v[0] + b*v0.v[2]; \ + p.v[1] = -a*v1.v[0] + b*v1.v[2]; \ + if (p.v[0] < p.v[1]) { min = p.v[0]; max = p.v[1]; } else { min = p.v[1]; max = p.v[0]; } \ + rad = fa * half_size -> v[0] + fb * half_size -> v[2]; \ + if (min > rad || max < -rad) return 0; + +/* ======================== Z-tests ======================== */ +/** + * + */ +#define AXISTEST_Z12(a, b, fa, fb) \ + p.v[1] = a*v1.v[0] - b*v1.v[1]; \ + p.v[2] = a*v2.v[0] - b*v2.v[1]; \ + if (p.v[2] < p.v[1]) { min = p.v[2]; max = p.v[1]; } else { min = p.v[1]; max = p.v[2]; } \ + rad = fa * half_size -> v[0] + fb * half_size -> v[1]; \ + if (min > rad || max < -rad) return 0; + +/** + * + */ +#define AXISTEST_Z0(a, b, fa, fb) \ + p.v[0] = a*v0.v[0] - b*v0.v[1]; \ + p.v[1] = a*v1.v[0] - b*v1.v[1]; \ + if (p.v[0] < p.v[1]) { min = p.v[0]; max = p.v[1]; } else { min = p.v[1]; max = p.v[0]; } \ + rad = fa * half_size -> v[0] + fb * half_size -> v[1]; \ + if (min > rad || max < -rad) return 0; + +#endif + +/** @} */ diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/include/kdtree.h blender/source/blender/render/intern/include/kdtree.h --- blender-cvs/source/blender/render/intern/include/kdtree.h 1969-12-31 19:00:00.000000000 -0500 +++ blender/source/blender/render/intern/include/kdtree.h 2006-06-14 03:35:34.000000000 -0400 @@ -0,0 +1,62 @@ +/* T I E . H + * + * @file kdtree.h + * + * BRL-CAD + * + * Copyright (c) 2002-2006 United States Government as represented by + * the U.S. Army Research Laboratory. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2 of + * the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this file; see the file named COPYING for more + * information. + * + * Comments - + * KD-Tree Routines + * + * Author - + * Justin L. Shumaker + * + * Source - + * The U. S. Army Research Laboratory + * Aberdeen Proving Ground, Maryland 21005-5068 USA + * + * $Id: kdtree.h,v 1.3 2006/01/18 06:46:11 brlcad Exp $ + */ +/** addtogroup libtie */ /** @{ */ + +#ifndef _KDTREE_H +#define _KDTREE_H + +#include "define.h" +#include "struct.h" + +#ifdef __cplusplus +extern "C" { +#endif + + +void tie_kdtree_free(tie_t *tie); +void tie_kdtree_cache_free(tie_t *tie, void **cache); +void tie_kdtree_cache_load(tie_t *tie, void *cache); +void tie_kdtree_prep(tie_t *tie); +extern tfloat TIE_PREC; + + +#ifdef __cplusplus +} +#endif + +#endif + +/** @} */ diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/include/render_types.h blender/source/blender/render/intern/include/render_types.h --- blender-cvs/source/blender/render/intern/include/render_types.h 2006-06-16 15:27:11.000000000 -0400 +++ blender/source/blender/render/intern/include/render_types.h 2006-06-14 13:01:47.000000000 -0400 @@ -40,6 +40,8 @@ #include "RE_pipeline.h" #include "RE_shader_ext.h" /* TexResult, ShadeResult, ShadeInput */ +#include "tie.h" /* triangle intersection engine */ + struct Object; struct MemArena; struct VertTableNode; @@ -144,6 +146,7 @@ /* octree tables and variables for raytrace */ Octree oc; + tie_t tie; /* use this instead of R.r.cfra */ float cfra; @@ -246,7 +249,7 @@ typedef struct VlakRen { - struct VertRen *v1, *v2, *v3, *v4; + struct VertRen *v1, *v2, *v3, *v4; /* 1 2 3 first tri, 1 3 4 second */ unsigned int lay; float n[3]; struct Material *mat; diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/include/struct.h blender/source/blender/render/intern/include/struct.h --- blender-cvs/source/blender/render/intern/include/struct.h 1969-12-31 19:00:00.000000000 -0500 +++ blender/source/blender/render/intern/include/struct.h 2006-06-14 03:35:34.000000000 -0400 @@ -0,0 +1,136 @@ +/* S T R U C T. H + * BRL-CAD + * + * Copyright (c) 2002-2006 United States Government as represented by + * the U.S. Army Research Laboratory. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2 of + * the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this file; see the file named COPYING for more + * information. + */ +/** @file struct.h + * S T R U C T . H + * + * Triangle Intersection Engine Structures + * + * Author - + * Justin L. Shumaker + * + * Source - + * The U. S. Army Research Laboratory + * Aberdeen Proving Ground, Maryland 21005-5068 USA + * + * $Id: struct.h,v 1.9 2006/01/18 06:46:11 brlcad Exp $ + */ +/** @addtogroup libtie */ /** @{ */ + +#ifndef _TIE_STRUCT_H +#define _TIE_STRUCT_H + +#include "define.h" +#include +#if TIE_RAY_ACCOUNTING +# include +#endif + +/** + * @struct TIE_3_s struct.h + * This is a 3-tuple + */ +typedef struct TIE_3_s { + tfloat v[3]; +} TIE_3; + +/** + * @struct tie_ray_s struct.h + * All the information about a ray. + * + * includes position, direction and depth + */ +typedef struct tie_ray_s { + TIE_3 pos; /* Position */ + TIE_3 dir; /* Direction */ + short depth; /* Depth */ + short kdtree_depth; +} tie_ray_t; + +/** + * @struct tie_id_s struct.h + * + * Ray intersection data (id) + * + * point, normal, distance, and barycentric coordinates + */ +typedef struct tie_id_s { + TIE_3 pos; /* Point */ + TIE_3 norm; /* Normal */ + tfloat dist; /* Distance */ + tfloat alpha; /* Barycentric Coordinate Alpha */ + tfloat beta; /* Barycentric Coordinate Beta */ +} tie_id_t; + +/** + * @struct tie_tri_s struct.h + * + * Everything you need to know about a triangle + */ +typedef struct tie_tri_s { + TIE_3 data[3]; /* 36-bytes, Data[0] = Point, Data[1] = Normal, Data[2] = DotProduct, VectorU, VectorV */ + tfloat *v; /* 8-bytes */ + void *ptr; /* 4-bytes */ +} tie_tri_t; + +/** + * @struct tie_kdtree_s struct.h + * + * The binary space partitioning tree + */ +typedef struct tie_kdtree_s { + tfloat axis; + void *data; +} tie_kdtree_t; + +/** + * + */ +typedef struct tie_geom_s { + tie_tri_t **tri_list; + int tri_num; +} tie_geom_t; + +/** + * + */ +typedef struct tie_stack_s { + tie_kdtree_t *node; + tfloat near; + tfloat far; +} tie_stack_t; + +/** + * + */ +typedef struct tie_s { + uint64_t rays_fired; + tie_kdtree_t *kdtree; + int max_depth; /* Maximum depth allowed for given geometry */ + TIE_3 min, max; + unsigned int tri_num; + unsigned int tri_num_alloc; + tie_tri_t *tri_list; + int stat; /* used for testing various statistics */ +} tie_t; + +#endif + +/** @} */ diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/include/tie.h blender/source/blender/render/intern/include/tie.h --- blender-cvs/source/blender/render/intern/include/tie.h 1969-12-31 19:00:00.000000000 -0500 +++ blender/source/blender/render/intern/include/tie.h 2006-06-14 03:35:34.000000000 -0400 @@ -0,0 +1,65 @@ +/* T I E . H + * + * @file tie.h + * + * BRL-CAD + * + * Copyright (c) 2002-2006 United States Government as represented by + * the U.S. Army Research Laboratory. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2 of + * the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this file; see the file named COPYING for more + * information. + * + * Comments - + * Triangle Intersection Engine Header + * + * Author - + * Justin L. Shumaker + * + * Source - + * The U. S. Army Research Laboratory + * Aberdeen Proving Ground, Maryland 21005-5068 USA + * + * $Id: tie.h,v 1.6 2006/01/18 06:46:11 brlcad Exp $ + */ +/** addtogroup libtie */ /** @{ */ + +#ifndef _TIE_H +#define _TIE_H + +#include "define.h" +#include "struct.h" +#include "kdtree.h" + +#ifdef __cplusplus +extern "C" { +#endif + + +void tie_init(tie_t *tie, unsigned int tri_num); +void tie_free(tie_t *tie); +void tie_prep(tie_t *tie); +void* tie_work(tie_t *tie, tie_ray_t *ray, tie_id_t *id, + void *(*hitfunc)(tie_ray_t*, tie_id_t*, tie_tri_t*, void *ptr), + void *ptr); +void tie_push(tie_t *tie, TIE_3 *tlist, int tnum, void *plist, int pstride); + + +#ifdef __cplusplus +} +#endif + +#endif + +/** @} */ diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/source/convertblender.c blender/source/blender/render/intern/source/convertblender.c --- blender-cvs/source/blender/render/intern/source/convertblender.c 2006-06-16 15:27:12.000000000 -0400 +++ blender/source/blender/render/intern/source/convertblender.c 2006-06-16 03:47:39.000000000 -0400 @@ -2933,7 +2933,8 @@ re->scene->world->aosphere= NULL; } - if(re->r.mode & R_RAYTRACE) freeoctree(re); + //if(re->r.mode & R_RAYTRACE) freeoctree(re); + if(re->r.mode & R_RAYTRACE) freekdtree(re); re->totvlak=re->totvert=re->totlamp=re->tothalo= 0; re->i.convertdone= 0; @@ -3352,7 +3353,8 @@ /* octree */ if(!re->test_break()) { if(re->r.mode & R_RAYTRACE) { - makeoctree(re); + /*makeoctree(re); */ + makekdtree(re); } } /* ENVIRONMENT MAPS */ diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/source/kdtree.c blender/source/blender/render/intern/source/kdtree.c --- blender-cvs/source/blender/render/intern/source/kdtree.c 1969-12-31 19:00:00.000000000 -0500 +++ blender/source/blender/render/intern/source/kdtree.c 2006-06-16 02:51:44.000000000 -0400 @@ -0,0 +1,929 @@ +/* K D T R E E . C +* +* @file tie.c +* +* BRL-CAD +* +* Copyright (c) 2002-2006 United States Government as represented by +* the U.S. Army Research Laboratory. +* +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public License +* as published by the Free Software Foundation; either version 2 of +* the License, or (at your option) any later version. +* +* This library is distributed in the hope that it will be useful, but +* WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +* Library General Public License for more details. +* +* You should have received a copy of the GNU Lesser General Public +* License along with this file; see the file named COPYING for more +* information. +* +* @brief support routines for shooting at triangles +* Comments - +* KD-Tree Routines +* +* Author - +* Justin L. Shumaker +* +* Source - +* The U. S. Army Research Laboratory +* Aberdeen Proving Ground, Maryland 21005-5068 USA +* +* $Id: kdtree.c,v 1.11 2006/01/18 06:46:11 brlcad Exp $ +*/ +/** @addtogroup libtie + * @{ + */ + +#include "kdtree.h" +#include +#include +#include +#include +#include + +#include "MEM_guardedalloc.h" + +tfloat TIE_PREC; + +/************************************************************* + **************** PRIVATE FUNCTIONS ************************** + *************************************************************/ + + +static void tie_kdtree_free_node(tie_kdtree_t *node) +{ + tie_kdtree_t *node_aligned = (tie_kdtree_t *)((intptr_t)node & ~0x7L); + + if (((intptr_t)(node_aligned->data)) & 0x4) { + /* Node Data is KDTREE Children, Recurse */ + tie_kdtree_free_node(&((tie_kdtree_t *)(((intptr_t)(node_aligned->data)) & ~0x7L))[0]); + tie_kdtree_free_node(&((tie_kdtree_t *)(((intptr_t)(node_aligned->data)) & ~0x7L))[1]); + MEM_freeN((void*)((intptr_t)(node_aligned->data) & ~0x7L)); + } else { + /* This node points to a geometry node, free it */ + MEM_freeN(((tie_geom_t *)((intptr_t)(node_aligned->data) & ~0x7L))->tri_list); + MEM_freeN((void *)((intptr_t)(node_aligned->data) & ~0x7L)); + } +} + + +static void tie_kdtree_cache_free_node(tie_t *tie, tie_kdtree_t *node, void **cache) +{ + tie_kdtree_t *node_aligned = (tie_kdtree_t *)((intptr_t)node & ~0x7L); + unsigned int size, mem, tri_num, i, tri_ind; + unsigned char type, split; + + memcpy(&size, *cache, sizeof(unsigned int)); + memcpy(&mem, &((char *)*cache)[sizeof(unsigned int)], sizeof(unsigned int)); + /* + * If the available size for this cache is under 1MB, then grow it 4MB larger. + * This reduces the number of realloc's required. Makes things much much faster + * on systems like FreeBSD that do a full copy for each realloc. + */ + if (mem - size < 1 << 20) { + mem += 1 << 23; + memcpy(&((char *)*cache)[sizeof(unsigned int)], &mem, sizeof(unsigned int)); + *cache = MEM_reallocN(*cache, mem); + } + + if (((intptr_t)(node_aligned->data)) & 0x4) { + /* Create a KD-Tree Node in the cache */ + type = 0; + memcpy(&((char *)*cache)[size], &type, 1); + size += 1; + memcpy(&((char *)*cache)[size], &(node_aligned->axis), sizeof(tfloat)); + size += sizeof(tfloat); + split = ((intptr_t)(node_aligned->data)) & 0x3; + memcpy(&((char *)*cache)[size], &split, 1); + size += 1; + + /* Update size of cache */ + memcpy(*cache, &size, sizeof(unsigned int)); + + /* Node Data is KDTREE Children, Recurse */ + tie_kdtree_cache_free_node(tie, &((tie_kdtree_t *)(((intptr_t)(node_aligned->data)) & ~0x7L))[0], cache); + tie_kdtree_cache_free_node(tie, &((tie_kdtree_t *)(((intptr_t)(node_aligned->data)) & ~0x7L))[1], cache); + MEM_freeN((void *)((intptr_t)(node_aligned->data) & ~0x7L)); + } else { + tri_num = ((tie_geom_t *)((intptr_t)(node_aligned->data) & ~0x7L))->tri_num; + type = 1; + memcpy(&((char *)*cache)[size], &type, 1); + size += 1; + memcpy(&((char *)*cache)[size], &tri_num, sizeof(unsigned int)); + + size += sizeof(unsigned int); + + for (i = 0; i < tri_num; i++) { + /* + * Pointer subtraction gives us the index of the triangle since the block of memory + * that the triangle exists in is contiguous memory. + */ + tri_ind = ((tie_geom_t *)((intptr_t)(node_aligned->data) & ~0x7L))->tri_list[i] - &tie->tri_list[0]; + memcpy(&((char *)*cache)[size], &tri_ind, sizeof(unsigned int)); + size += sizeof(unsigned int); + } + + /* Update size of cache */ + memcpy(*cache, &size, sizeof(unsigned int)); + + /* This node points to a geometry node, free it */ + MEM_freeN(((tie_geom_t *)((intptr_t)(node_aligned->data) & ~0x7L))->tri_list); + MEM_freeN((void *)((intptr_t)(node_aligned->data) & ~0x7L)); + } +} + + +static void tie_kdtree_prep_head(tie_t *tie, tie_tri_t *tri_list, int tri_num) +{ + tie_geom_t *g; + TIE_3 min, max; + int i; + + + if (!tri_num) + return ; + + /* Insert all triangles into the Head Node */ + if (!tie->kdtree) { + tie->kdtree = (tie_kdtree_t *)MEM_mallocN(sizeof(tie_kdtree_t), "kdtree"); + tie->kdtree->data = (void *)MEM_mallocN(sizeof(tie_geom_t), "kdtree data"); + g = ((tie_geom_t *)(tie->kdtree->data)); + g->tri_num = 0; + + math_bbox(tie->min, tie->max, tri_list[0].data[0], tri_list[0].data[1], tri_list[0].data[2]); + + g->tri_list = (tie_tri_t **)MEM_mallocN(sizeof(tie_tri_t *) * tri_num, "kdtree tri_list"); + + /* form bounding box of scene */ + for (i = 0; i < tri_num; i++) { + g->tri_list[i] = &tri_list[i]; + + /* Get Bounding Box of Triangle */ + math_bbox(min, max, tri_list[i].data[0], tri_list[i].data[1], tri_list[i].data[2]); + + /* Check to see if defines a new Max or Min point */ + math_vec_min(tie->min, min); + math_vec_max(tie->max, max); + /* printf("Box: [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n", tie->min.v[0], tie->min.v[1], tie->min.v[2], tie->max.v[0], tie->max.v[1], tie->max.v[2]); */ + } + + ((tie_geom_t *)(tie->kdtree->data))->tri_num = tri_num; + } +} + + +static int tie_kdtree_tri_box_overlap(TIE_3 *center, TIE_3 *half_size, TIE_3 triverts[3]) +{ + /* + * use separating axis theorem to test overlap between triangle and box + * need to test for overlap in these directions: + * 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle + * we do not even need to test these) + * 2) normal of the triangle + * 3) crossproduct(edge from tri, {x,y,z}-directin) + * this gives 3x3=9 more tests + */ + TIE_3 v0, v1, v2, normal, e0, e1, e2, fe, p; + tfloat min, max, d, t, rad; + + /* move everything so that the boxcenter is in (0,0,0) */ + math_vec_sub(v0, triverts[0], (*center)); + math_vec_sub(v1, triverts[1], (*center)); + math_vec_sub(v2, triverts[2], (*center)); + + /* + * First test overlap in the {x,y,z}-directions + * find min, max of the triangle each direction, and test for overlap in + * that direction -- this is equivalent to testing a minimal AABB around + * the triangle against the AABB + */ + + /* test in X-direction */ + math_min3(min, v0.v[0], v1.v[0], v2.v[0]); + math_max3(max, v0.v[0], v1.v[0], v2.v[0]); + if (min > half_size->v[0] || max < -half_size->v[0]) + return 0; + + /* test in Y-direction */ + math_min3(min, v0.v[1], v1.v[1], v2.v[1]); + math_max3(max, v0.v[1], v1.v[1], v2.v[1]); + if (min > half_size->v[1] || max < -half_size->v[1]) + return 0; + + /* test in Z-direction */ + math_min3(min, v0.v[2], v1.v[2], v2.v[2]); + math_max3(max, v0.v[2], v1.v[2], v2.v[2]); + if (min > half_size->v[2] || max < -half_size->v[2]) + return 0; + + /* compute triangle edges */ + math_vec_sub(e0, v1, v0); /* tri edge 0 */ + math_vec_sub(e1, v2, v1); /* tri edge 1 */ + math_vec_sub(e2, v0, v2); /* tri edge 2 */ + + /* Perform the 9 tests */ + fe.v[0] = fabs(e0.v[0]); + fe.v[1] = fabs(e0.v[1]); + fe.v[2] = fabs(e0.v[2]); + + AXISTEST_X01(e0.v[2], e0.v[1], fe.v[2], fe.v[1]); + AXISTEST_Y02(e0.v[2], e0.v[0], fe.v[2], fe.v[0]); + AXISTEST_Z12(e0.v[1], e0.v[0], fe.v[1], fe.v[0]); + + fe.v[0] = fabs(e1.v[0]); + fe.v[1] = fabs(e1.v[1]); + fe.v[2] = fabs(e1.v[2]); + + AXISTEST_X01(e1.v[2], e1.v[1], fe.v[2], fe.v[1]); + AXISTEST_Y02(e1.v[2], e1.v[0], fe.v[2], fe.v[0]); + AXISTEST_Z0(e1.v[1], e1.v[0], fe.v[1], fe.v[0]); + + fe.v[0] = fabs(e2.v[0]); + fe.v[1] = fabs(e2.v[1]); + fe.v[2] = fabs(e2.v[2]); + + AXISTEST_X2(e2.v[2], e2.v[1], fe.v[2], fe.v[1]); + AXISTEST_Y1(e2.v[2], e2.v[0], fe.v[2], fe.v[0]); + AXISTEST_Z12(e2.v[1], e2.v[0], fe.v[1], fe.v[0]); + + /* + * Test if the box intersects the plane of the triangle + * compute plane equation of triangle: normal*x+d=0 + */ + math_vec_cross(normal, e0, e1); + math_vec_dot(d, normal, v0); /* plane eq: normal . x + d = 0 */ + + p.v[0] = normal.v[0] > 0 ? half_size->v[0] : -half_size->v[0]; + p.v[1] = normal.v[1] > 0 ? half_size->v[1] : -half_size->v[1]; + p.v[2] = normal.v[2] > 0 ? half_size->v[2] : -half_size->v[2]; + + math_vec_dot(t, normal, p); + return t >= d ? 1 : 0; +} + + +static void tie_kdtree_build(tie_t *tie, tie_kdtree_t *node, int depth, TIE_3 min, TIE_3 max, int node_a, int node_b) +{ + tie_geom_t *child[2], *node_gd = (tie_geom_t *)(node->data); + TIE_3 cmin[2], cmax[2], center[2], half_size[2]; + int i, j, n, split, cnt[2]; + +#if 0 + /* if(depth >= 26) */ + printf("%f %f %f %f %f %f\n", min.v[0], min.v[1], min.v[2], max.v[0], max.v[1], max.v[2]); +#endif + + /* Terminating criteria for KDTREE subdivision */ + if (node_gd->tri_num <= TIE_KDTREE_NODE_MAX || depth > tie->max_depth) { + /* tie->stat++; */ + tie->stat += node_gd->tri_num; +#if 0 + + if (node_gd->tri_num > tie->stat) + tie->stat = node_gd->tri_num; + if (node_gd->tri_num > tie->stat) { + tie->stat = node_gd->tri_num; + printf("depth: %d, tris: %d\n", depth, node_gd->tri_num); + printf("%f %f %f %f %f %f\n", min.v[0], min.v[1], min.v[2], max.v[0], max.v[1], max.v[2]); + } + exit(0); +#endif + + return ; + } + +#if 0 + { + /********************** + * MID-SPLIT ALGORITHM * + ***********************/ + TIE_3 vec; + + /* Left Child */ + cmin[0] = min; + cmax[0] = max; + + /* Right Child */ + cmin[1] = min; + cmax[1] = max; + + math_vec_add(center[0], max, min); + math_vec_mul_scalar(center[0], center[0], 0.5); + + /* Split along largest Axis to keep node sizes relatively cube-like (Naive) */ + math_vec_sub(vec, max, min); + + /* Determine the largest Axis */ + if (vec.v[0] >= vec.v[1] && vec.v[0] >= vec.v[2]) + { + cmax[0].v[0] = center[0].v[0]; + cmin[1].v[0] = center[0].v[0]; + node->axis = center[0].v[0]; + split = 0; + } else if (vec.v[1] >= vec.v[0] && vec.v[1] >= vec.v[2]) + { + cmax[0].v[1] = center[0].v[1]; + cmin[1].v[1] = center[0].v[1]; + node->axis = center[0].v[1]; + split = 1; + } else + { + cmax[0].v[2] = center[0].v[2]; + cmin[1].v[2] = center[0].v[2]; + node->axis = center[0].v[2]; + split = 2; + } + } +#else + { + /**************************************** + * Justin's Aggressive KD-Tree Algorithm * + *****************************************/ + int slice[3][MAX_SLICES + MIN_SLICES], gap[3][2], active, split_slice; + int side[3][MAX_SLICES + MIN_SLICES][2], d, s, k, smax[3], smin, slice_num; + tfloat coef[3][MAX_SLICES + MIN_SLICES], split_coef, beg, end, d_min, d_max; + tie_tri_t *tri; + + /* + * Calculate number of slices to use as a function of triangle density. + * Slices as a function of relative node size does not work so well. + */ + slice_num = MIN_SLICES + MAX_SLICES * ((tfloat)node_gd->tri_num / (tfloat)tie->tri_num); + + for (d = 0; d < 3; d++) + { + /* + * Optimization: Walk each triangle and find the min and max for the given dimension + * of the complete triangle list. This will tell us what slices we needn't bother + * doing any computations for. + */ + for (i = 0; i < node_gd->tri_num; i++) { + tri = node_gd->tri_list[i]; + /* Set min anx max */ + math_min3(tri->v[0], tri->data[0].v[d], tri->data[1].v[d], tri->data[2].v[d]); + math_max3(tri->v[1], tri->data[0].v[d], tri->data[1].v[d], tri->data[2].v[d]); + + /* Clamp to node AABB */ + if (tri->v[0] < min.v[d]) + tri->v[0] = min.v[d]; + if (tri->v[1] > max.v[d]) + tri->v[1] = max.v[d]; + + if (i == 0 || tri->v[0] < d_min) + d_min = tri->v[0]; + + if (i == 0 || tri->v[1] > d_max) + d_max = tri->v[1]; + } + + for (k = 0; k < slice_num; k++) { + slice[d][k] = 0; + side[d][k][0] = 0; + side[d][k][1] = 0; + + /* Left Child */ + cmin[0] = min; + cmax[0] = max; + + /* Right Child */ + cmin[1] = min; + cmax[1] = max; + + /* construct slices so as not to use the boundaries as slices */ + coef[d][k] = ((tfloat)k / (tfloat)(slice_num - 1)) * (tfloat)(slice_num - 2) / (tfloat)slice_num + (tfloat)1 / (tfloat)slice_num; + cmax[0].v[d] = min.v[d] * (1.0 - coef[d][k]) + max.v[d] * coef[d][k]; + cmin[1].v[d] = cmax[0].v[d]; + + if (cmax[0].v[d] < d_min || cmax[0].v[d] > d_max) + continue; + + for (n = 0; n < 2; n++) { + math_vec_add(center[n], cmax[n], cmin[n]); + math_vec_mul_scalar(center[n], center[n], 0.5); + math_vec_sub(half_size[n], cmax[n], cmin[n]); + math_vec_mul_scalar(half_size[n], half_size[n], 0.5); + } + + for (i = 0; i < node_gd->tri_num; i++) { + /* + * Optimization: If the points for the triangle of the dimension being tested + * do not span the cutting plane, then do not bother with the next test. + */ + if ((node_gd->tri_list[i]->data[0].v[d] > cmax[0].v[d] && + node_gd->tri_list[i]->data[1].v[d] > cmax[0].v[d] && + node_gd->tri_list[i]->data[2].v[d] > cmax[0].v[d]) || + (node_gd->tri_list[i]->data[0].v[d] < cmax[0].v[d] && + node_gd->tri_list[i]->data[1].v[d] < cmax[0].v[d] && + node_gd->tri_list[i]->data[2].v[d] < cmax[0].v[d])) + continue; + + /* Check that the triangle is in both node A and B for it to span. */ + s = 0; + for (n = 0; n < 2; n++) { + /* + * Check to see if any triangle points are inside of the node before + * spending alot of cycles on the full blown triangle box overlap + */ + for (j = 0; j < 3; j++) + if (node_gd->tri_list[i]->data[j].v[0] > cmin[n].v[0] && + node_gd->tri_list[i]->data[j].v[0] < cmax[n].v[0] && + node_gd->tri_list[i]->data[j].v[1] > cmin[n].v[1] && + node_gd->tri_list[i]->data[j].v[1] < cmax[n].v[1] && + node_gd->tri_list[i]->data[j].v[2] > cmin[n].v[2] && + node_gd->tri_list[i]->data[j].v[2] < cmax[n].v[2]) { + j = 4; + } + + if (j == 5) { + s++; + side[d][k][n]++; + } else { + if (tie_kdtree_tri_box_overlap(¢er[n], &half_size[n], node_gd->tri_list[i]->data)) { + s++; + side[d][k][n]++; + } + } + } + + if (s == 2) + slice[d][k]++; + } + } + } + + /* Store the max value from each of the 3 Slice arrays */ + for (d = 0; d < 3; d++) + { + smax[d] = 0; + for (k = 0; k < slice_num; k++) { + if (slice[d][k] > smax[d]) + smax[d] = slice[d][k]; + } + } + + /* + * To eliminate "empty" areas, build a list of spans where geometric complexity is + * less than MIN_SPAN of the overal nodes size and then selecting the splitting plane + * the corresponds to the span slice domain nearest the center to bias towards a balanced tree + */ + + for (d = 0; d < 3; d++) + { + gap[d][0] = 0; + gap[d][1] = 0; + beg = 0; + end = 0; + active = 0; + + for (k = 0; k < slice_num; k++) { + /* printf("slice[%d][%d]: %d < %d\n", d, k, slice[d][k], (int)(MIN_DENSITY * (tfloat)smax[d])); */ + if (slice[d][k] < (int)(MIN_DENSITY * (tfloat)smax[d])) { + if (!active) { + active = 1; + beg = k; + end = k; + } else { + end = k; + } + } else { + if (active) { + if (end - beg > gap[d][1] - gap[d][0]) { + gap[d][0] = beg; + gap[d][1] = end; + } + } + active = 0; + } + } + + if (active) + if (end - beg > gap[d][1] - gap[d][0]) { + gap[d][0] = beg; + gap[d][1] = end; + } + } + +#if 0 + printf("gap_x: %d->%d = %d\n", gap[0][0], gap[0][1], gap[0][1] - gap[0][0]); + printf("gap_y: %d->%d = %d\n", gap[1][0], gap[1][1], gap[1][1] - gap[1][0]); + printf("gap_z: %d->%d = %d\n", gap[2][0], gap[2][1], gap[2][1] - gap[2][0]); +#endif + + /* + * If there is a gap atleast MIN_SPAN in side wrt the nodes dimension size + * then use the nearest edge of the gap to 0.5 as the splitting plane, + * Use the the gap with the largest span. + * If no gaps are found meeting the criteria then weight the span values to + * bias towards a balanced kd-tree and choose the minima of that weighted curve. + */ + + /* Largest gap */ + d = 0; + if (gap[1][1] - gap[1][0] > gap[d][1] - gap[d][0]) + d = 1; + if (gap[2][1] - gap[2][0] > gap[d][1] - gap[d][0]) + d = 2; + + /* + * Largest gap found must meet MIN_SPAN requirements + * There must be atleast 500 triangles or we don't bother. + * Lower triangle numbers means there is a higher probability that + * triangles lack any sort of coherent structure. + */ + if ((tfloat)(gap[d][1] - gap[d][0]) / (tfloat)slice_num > MIN_SPAN && node_gd->tri_num > 500) + { + /* printf("choosing slice[%d]: %d->%d :: %d tris\n", d, gap[d][0], gap[d][1], node_gd->tri_num); */ + split = d; + if (abs(gap[d][0] - slice_num / 2) < abs(gap[d][1] - slice_num / 2)) { + /* choose gap[d][0] as splitting plane */ + split_coef = ((tfloat)gap[d][0] / (tfloat)(slice_num - 1)) * (tfloat)(slice_num - 2) / (tfloat)slice_num + (tfloat)1 / (tfloat)slice_num; + split_slice = gap[d][0]; + } else { + /* choose gap[d][1] as splitting plane */ + split_coef = ((tfloat)gap[d][1] / (tfloat)(slice_num - 1)) * (tfloat)(slice_num - 2) / (tfloat)slice_num + (tfloat)1 / (tfloat)slice_num; + split_slice = gap[d][1]; + } + } else + { + /* + * Weight the slices based on a heuristic driven linear scaling function to bias values + * towards the center as more desirable. This solves the case of a partially linear graph + * to prevent marching to determine a desirable splitting point. If this section of code + * is being executed it's typically because most 'empty space' has now been eliminated + * and/or the resulting geometry is now losing structure as the smaller cells are being + * created, i.e dividing a fraction of a wing-nut instead of an engine-block. + */ + for (d = 0; d < 3; d++) { + for (k = 0; k < slice_num; k++) { + slice[d][k] += fabs(coef[d][k] - 0.5) * SCALE_COEF * smax[d]; + /* printf("%.3f %d\n", coef[d][k], slice[d][k]); */ + } + } + + /* Choose the slice with the graphs minima as the splitting plane. */ + split = 0; + smin = tie->tri_num; + split_coef = 0.5; + for (d = 0; d < 3; d++) { + for (k = 0; k < slice_num; k++) { + if (slice[d][k] < smin) { + split_coef = coef[d][k]; + split = d; + split_slice = k; + smin = slice[d][k]; + } + } + } + + /* + * If the dimension chosen to split along has a value of 0 for the maximum value + * then the geometry was aligned such that it fell undetectable between the slices + * and therefore was not picked up by the marching slices. In the event that this + * happens, choose to naively split along the middle as this last ditch decision + * will give better results than the algorithm naively picking the first of the + * the slices forming these irregular, short followed by a long box, splits. + */ + if (smax[split] == 0) { + split_coef = coef[split][slice_num / 2]; + split_slice = slice_num / 2; + } + } + + /* + * Lastly, after we have supposedly chosen the most ideal splitting point, + * check to see that the subdivision that is about to take place is worth + * doing. In other words, if both children have the same number of triangles + * as the parent does then stop. + */ + if (side[split][split_slice][0] == node_gd->tri_num && side[split][split_slice][1] == node_gd->tri_num) + { + tie->stat += node_gd->tri_num; + return ; + } + +#if 0 + if (side[split][split_slice][0] == node_a && side[split][split_slice][1] == node_b) + { + if (node_gd->tri_num < 10) + return ; + /* printf("%f %f %f %f %f %f\n", min.v[0], min.v[1], min.v[2], max.v[0], max.v[1], max.v[2]); */ + /* printf("moo: %d - %d\n", depth, node_gd->tri_num); */ + } +#endif + + +#if 0 + printf("winner: depth: %d, dim = %d, smin = %d, coef: %.3f\n", depth, split, smin, split_coef); + printf("winner: min: %.3f %.3f %.3f, max: %.3f %.3f %.3f, tris: %d\n", min.v[0], min.v[1], min.v[2], max.v[0], max.v[1], max.v[2], node_gd->tri_num); +#endif + + /* Based on the winner, construct the two child nodes */ + /* Left Child */ + cmin[0] = min; + cmax[0] = max; + + /* Right Child */ + cmin[1] = min; + cmax[1] = max; + + cmax[0].v[split] = min.v[split] * (1.0 - split_coef) + max.v[split] * split_coef; + cmin[1].v[split] = cmax[0].v[split]; + node->axis = cmax[0].v[split]; + } +#endif + + + /* Allocate 2 children nodes for the parent node */ + node->data = (void *)MEM_mallocN(2 * sizeof(tie_kdtree_t), "kdtree build"); + + ((tie_kdtree_t *)(node->data))[0].data = MEM_mallocN(sizeof(tie_geom_t), "kdtree node data"); + ((tie_kdtree_t *)(node->data))[1].data = MEM_mallocN(sizeof(tie_geom_t), "kdtree node data"); + + /* Initialize Triangle List */ + child[0] = ((tie_geom_t *)(((tie_kdtree_t *)(node->data))[0].data)); + child[1] = ((tie_geom_t *)(((tie_kdtree_t *)(node->data))[1].data)); + + child[0]->tri_list = (tie_tri_t **)MEM_mallocN(sizeof(tie_tri_t *) * node_gd->tri_num, "kdtree tri_list"); + child[0]->tri_num = 0; + + child[1]->tri_list = (tie_tri_t **)MEM_mallocN(sizeof(tie_tri_t *) * node_gd->tri_num, "kdtree tri_list"); + child[1]->tri_num = 0; + + + /* + * Determine if the triangles touch either of the two children nodes, + * if it does insert it into them respectively. + */ + for (n = 0; n < 2; n++) { + cnt[n] = 0; + + math_vec_add(center[n], cmax[n], cmin[n]); + math_vec_mul_scalar(center[n], center[n], 0.5); + math_vec_sub(half_size[n], cmax[n], cmin[n]); + math_vec_mul_scalar(half_size[n], half_size[n], 0.5); + + for (i = 0; i < node_gd->tri_num; i++) { + /* + * Check to see if any triangle points are inside of the node before + * spending alot of cycles on the full blown triangle box overlap + */ + for (j = 0; j < 3; j++) + if (node_gd->tri_list[i]->data[j].v[0] > cmin[n].v[0] && + node_gd->tri_list[i]->data[j].v[0] < cmax[n].v[0] && + node_gd->tri_list[i]->data[j].v[1] > cmin[n].v[1] && + node_gd->tri_list[i]->data[j].v[1] < cmax[n].v[1] && + node_gd->tri_list[i]->data[j].v[2] > cmin[n].v[2] && + node_gd->tri_list[i]->data[j].v[2] < cmax[n].v[2]) { + j = 4; + } + + if (j == 5) { + child[n]->tri_list[child[n]->tri_num++] = node_gd->tri_list[i]; + cnt[n]++; + } else { + if (tie_kdtree_tri_box_overlap(¢er[n], &half_size[n], node_gd->tri_list[i]->data)) { + child[n]->tri_list[child[n]->tri_num++] = node_gd->tri_list[i]; + cnt[n]++; + } + } + } + + /* Resize Tri List to actual ammount of memory used */ + child[n]->tri_list = (tie_tri_t **)MEM_reallocN(child[n]->tri_list, sizeof(tie_tri_t *) * child[n]->tri_num); + } + + /* + * Now that the triangles have been propogated to the appropriate child nodes, + * free the triangle list on this node. + */ + node_gd->tri_num = 0; + MEM_freeN(node_gd->tri_list); + MEM_freeN(node_gd); + + /* Push each child through the same process. */ + tie_kdtree_build(tie, &((tie_kdtree_t *)(node->data))[0], depth + 1, cmin[0], cmax[0], cnt[0], cnt[1]); + tie_kdtree_build(tie, &((tie_kdtree_t *)(node->data))[1], depth + 1, cmin[1], cmax[1], cnt[0], cnt[1]); + + /* Assign the splitting dimension to the node */ + /* If we've come this far then YES, this node DOES have child nodes, MARK it as so. */ + node->data = (void *)((intptr_t)(node->data) + split + 4); +} + + +/************************************************************* + **************** EXPORTED FUNCTIONS ************************* + *************************************************************/ + +/** + * Free up all the stuff associated with the kdtree + * + * All of the KDTREE nodes and triangles that we have allocated need to + * be freed in a controlled manner. This routine does that. + * + * @param tie pointer to a struct tie_t + * @return void + */ +void tie_kdtree_free(tie_t *tie) +{ + /* Free KDTREE Nodes */ + /* prevent tie from crashing when a tie_free() is called right after a tie_init() */ + if (tie->kdtree) + tie_kdtree_free_node(tie->kdtree); + MEM_freeN(tie->kdtree); +} + + +/** + * Free up all the stuff associated with the kdtree and build a + * cache as the data is freed. Building the cache while the data + * is freed allows the peak memory not to go any higher than it + * already is. If there were seprarate cache and free functions + * then the cache would exist in memory while the triangles and + * kd-tree were in memory thus severly limiting optimal memory + * usage. + * + * All of the KDTREE nodes and triangles that we have allocated need to + * be freed in a controlled manner. This routine does that. + * + * @param tie pointer to a struct tie_t and void **cache to store data + * @return void + */ +void tie_kdtree_cache_free(tie_t *tie, void **cache) +{ + unsigned int size; + + /* + * Free KDTREE Node + * Prevent tie from crashing when a tie_free() is called right after a tie_init() + */ + if (tie->kdtree) { + *cache = MEM_mallocN(2 * sizeof(unsigned int), "kdtree cache"); + size = 2 * sizeof(unsigned int); + memcpy(*cache, &size, sizeof(unsigned int)); + memcpy(&((char *)*cache)[sizeof(unsigned int)], &size, sizeof(unsigned int)); + + /* Build the cache */ + tie_kdtree_cache_free_node(tie, tie->kdtree, cache); + + /* Resize the array back to it's real value */ + memcpy(&size, *cache, sizeof(unsigned int)); + *cache = MEM_reallocN(*cache, size); + } + MEM_freeN(tie->kdtree); +} + + +void tie_kdtree_cache_load(tie_t *tie, void *cache) +{ + tie_kdtree_t *node, *temp_node, *stack[64]; + tie_geom_t *geom; + TIE_3 min, max; + unsigned int i, size, index, tri_ind, stack_ind; + char type, split; + + + if (!cache) + return ; + + memcpy(&size, cache, sizeof(unsigned int)); + /* Advance past the first (2) unsigned ints to the actualy data */ + index = 2 * sizeof(unsigned int); + stack_ind = 0; + + while (index < size) { + memcpy(&type, &((char *)cache)[index], 1); + index += 1; + + if (type) { + /* Geometry Node - Allocate a tie_geom_t and assign to node->data. */ + node->data = MEM_mallocN(sizeof(tie_geom_t), "kdtree node data"); + geom = (tie_geom_t *)node->data; + + memcpy(&(geom->tri_num), &((char *)cache)[index], sizeof(unsigned int)); + index += sizeof(unsigned int); + + geom->tri_list = (tie_tri_t **)MEM_mallocN(geom->tri_num * sizeof(tie_tri_t *), "kdtree tri_list"); + for (i = 0; i < geom->tri_num; i++) { + memcpy(&tri_ind, &((char *)cache)[index], sizeof(unsigned int)); + index += sizeof(unsigned int); + + /* Translate the numerical index to a pointer index into tie->tri_list. */ + geom->tri_list[i] = &tie->tri_list[0] + tri_ind; + } + + if (stack_ind) { + stack_ind--; + node = stack[stack_ind]; + } + } else { + /* KD-Tree Node */ + if (!tie->kdtree) { + tie->kdtree = (tie_kdtree_t *)MEM_mallocN(sizeof(tie_kdtree_t), "kdtree cache load"); + node = tie->kdtree; + } + + /* Assign splitting axis value */ + memcpy(&node->axis, &((char *)cache)[index], sizeof(tfloat)); + index += sizeof(tfloat); + + /* Get splitting plane */ + memcpy(&split, &((char *)cache)[index], 1); + index += 1; + + /* Allocate memory for 2 child nodes */ + node->data = MEM_mallocN(2 * sizeof(tie_kdtree_t), "kdtree cache load"); + + /* Push B on the stack and Process A */ + stack[stack_ind] = &((tie_kdtree_t *)node->data)[1]; + stack_ind++; + + /* Set the new current node */ + temp_node = node; + node = &((tie_kdtree_t *)node->data)[0]; + + /* + * Mask the splitting plane and mark it as a kdtree node + * using the lower bits of the ptr. + */ + temp_node->data = (void *)((intptr_t)(temp_node->data) + split + 4); + } + } + + /* form bounding box of scene */ + math_bbox(tie->min, tie->max, tie->tri_list[0].data[0], tie->tri_list[0].data[1], tie->tri_list[0].data[2]); + for (i = 0; i < tie->tri_num; i++) { + /* Get Bounding Box of Triangle */ + math_bbox(min, max, tie->tri_list[i].data[0], tie->tri_list[i].data[1], tie->tri_list[i].data[2]); + + /* Check to see if defines a new Max or Min point */ + math_vec_min(tie->min, min); + math_vec_max(tie->max, max); + } +} + + +/** + * Get ready to shoot rays at triangles + * + * Build the KDTREE tree for the triangles we have + * + * @param tie pointer to a struct tie_t which now has all the triangles in it + * @return void + */ +void tie_kdtree_prep(tie_t *tie) +{ + TIE_3 delta; + int already_built; + + + already_built = tie->kdtree ? 1 : 0; + + /* Set bounding volume and make head node a geometry node */ + if (!already_built) + tie_kdtree_prep_head(tie, tie->tri_list, tie->tri_num); + + if (!tie->kdtree) + return ; + + /* Trim KDTREE to number of actual triangles if it's not that size already. */ + if (!already_built) + ((tie_geom_t *)(tie->kdtree->data))->tri_list = (tie_tri_t **)MEM_reallocN(((tie_geom_t *)(tie->kdtree->data))->tri_list, sizeof(tie_tri_t *) * ((tie_geom_t *)(tie->kdtree->data))->tri_num); + + /* + * Compute Floating Fuzz Precision Value + * For now, take largest dimension as basis for TIE_PREC + */ + math_vec_sub(delta, tie->max, tie->min); + math_max3(TIE_PREC, delta.v[0], delta.v[1], delta.v[2]); +#if TIE_SINGLE_PREC + + TIE_PREC *= 0.000001; +#else + + TIE_PREC *= 0.000000000001; +#endif + + /* Grow the head node to avoid floating point fuzz in the building process with edges */ + math_vec_mul_scalar(delta, delta, 1.0); /* XXX */ + math_vec_sub(tie->min, tie->min, delta); + math_vec_add(tie->max, tie->max, delta); + + /* Compute Max Depth to allow the KD-Tree to grow to */ + tie->max_depth = (int)(TIE_KDTREE_DEPTH_K1 * (log(tie->tri_num) / log(2)) + TIE_KDTREE_DEPTH_K2); + printf("max_depth: %d\n", tie->max_depth); + + /* Build the KDTREE */ + if (!already_built) + tie_kdtree_build(tie, tie->kdtree, 0, tie->min, tie->max, 0, 0); + + printf("stat: %d\n", tie->stat); + tie->stat = 0; + + /* exit(0); */ /* uncomment to profile prep phase only */ +} + +/** @} */ diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/source/pipeline.c blender/source/blender/render/intern/source/pipeline.c --- blender-cvs/source/blender/render/intern/source/pipeline.c 2006-06-16 14:41:18.000000000 -0400 +++ blender/source/blender/render/intern/source/pipeline.c 2006-06-16 17:43:20.000000000 -0400 @@ -724,6 +724,7 @@ /* init some variables */ re->ycor= 1.0f; + re->tie.kdtree = NULL; return re; } diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/source/ray.c blender/source/blender/render/intern/source/ray.c --- blender-cvs/source/blender/render/intern/source/ray.c 2006-06-16 15:27:12.000000000 -0400 +++ blender/source/blender/render/intern/source/ray.c 2006-06-16 17:53:35.000000000 -0400 @@ -653,6 +653,107 @@ re->stats_draw(&re->i); } +void freekdtree(Render *re) +{ + + printf("Freeing kdtree\n"); + if (re->tie.kdtree) { + printf(" actually freeing\n"); + tie_free(&re->tie); + re->tie.kdtree = NULL; /* HACK HACK BUG SCREAM */ + } + fflush(stdout); +} + +/* The new kdtree method; first implementation is Twingy's code from BRL-CAD. + * Later I'll try my own tree which stores hyperrectangles rather than + * splitting planes and is a bit more space efficient. */ +void makekdtree(Render *re) +{ + tie_t *tie = &re->tie; + VlakRen *vlr = NULL; + //double lasttime= PIL_check_seconds_timer(); + int triangle_count = 0; + int v; + + printf("Making kdtree\n"); + if (re->tie.kdtree) { + printf(" but it already exists. bailing for now rather than reconstructing\n"); + fflush(stdout); + return; + } + fflush(stdout); + + /* Update the status bar */ + re->i.infostr = "Counting triangles"; + re->stats_draw(&re->i); + + /* Only for debug info */ + raycount = 0; + accepted = 0; + rejected = 0; + coherent_ray = 0; + + /* Count how many triangles we have */ + triangle_count = 0; + for(v=0;vtotvlak;v++) { + if((v & 255)==0) + vlr = re->blovl[v>>8]; + else + vlr++; + if((vlr->mat->mode & MA_TRACEBLE) && + ((vlr->mat->mode & MA_WIRE)==0)) { + triangle_count++; + if(vlr->v4) { + triangle_count++; + } + } + } + + re->i.infostr = "Initializing the kdtree"; + re->stats_draw(&re->i); + + /* Allocate enough space in the kdtree */ + tie_init(tie, triangle_count); + + /* Push each triangle into the triangle intersection engine */ + for(v=0;vtotvlak;v++) { + if((v & 255)==0) + vlr= re->blovl[v>>8]; + else + vlr++; + if((vlr->mat->mode & MA_TRACEBLE) && + ((vlr->mat->mode & MA_WIRE)==0)) { + float tlist[9]; + tlist[0] = vlr->v1->co[0]; + tlist[1] = vlr->v1->co[1]; + tlist[2] = vlr->v1->co[2]; + tlist[3] = vlr->v2->co[0]; + tlist[4] = vlr->v2->co[1]; + tlist[5] = vlr->v2->co[2]; + tlist[6] = vlr->v3->co[0]; + tlist[7] = vlr->v3->co[1]; + tlist[8] = vlr->v3->co[2]; + tie_push(tie, (TIE_3*)tlist, 1, vlr, 4); + if(vlr->v4) { + tlist[6] = vlr->v4->co[0]; + tlist[7] = vlr->v4->co[1]; + tlist[8] = vlr->v4->co[2]; + tie_push(tie, (TIE_3*)tlist, 1, vlr, 4); + } + } + } + + re->i.infostr = "Building the kdtree"; + re->stats_draw(&re->i); + + /* Now do the hard work an build the kdtree */ + tie_prep(tie); + + re->i.infostr = NULL; + re->stats_draw(&re->i); +} + /* ************ raytracer **************** */ /* only for self-intersecting test with current render face (where ray left) */ @@ -949,7 +1050,7 @@ VlakRen *vlr; short nr=0; OcVal *ov; - + if(is->mode==DDA_SHADOW) { vlr= no->v[0]; @@ -1142,7 +1243,8 @@ /* return 1: found valid intersection */ /* starts with is->vlrorig */ -static int d3dda(Isect *is) +/* replaced by d3dda kdtree traversal below */ +static int d3dda_old(Isect *is) { Node *no; OcVal ocval; @@ -1387,6 +1489,72 @@ return 0; } +/* callback for triangle intersection; invoked by tie_work when an intersection + * is found. */ +void *tie_hitfunc(tie_ray_t* ray, tie_id_t* id, tie_tri_t* tri, void *ptr) +{ + Isect *is = ptr; + + is->start[0] = ray->pos.v[0]; + is->start[1] = ray->pos.v[1]; + is->start[2] = ray->pos.v[2]; + + /* FIXME is->vlrisect? */ + /* FIXME is->* i.e. other members? */ + + /* we don't want to continue intersecting... other code will do that */ + /* clearly we should be smarter about this but for now this is just a + * proof of concept */ + return (void*) 1; +} + +/* return 1: found valid intersection */ +/* starts with is->vlrorig */ +/* now using the kdtree */ +/* mode DDA_MIRROR: get closest intersection */ +/* mode DDA_SHADOW: get any intersection */ +/* mode DDA_SHADOW_TRA: get any intersection */ +static int d3dda(Isect *is) +{ + /* init tie_ray struct */ + tie_ray_t ray; + tie_id_t id; + ray.pos.v[0] = is->start[0]; + ray.pos.v[1] = is->start[1]; + ray.pos.v[2] = is->start[2]; + ray.dir.v[0] = is->vec[0]; + ray.dir.v[1] = is->vec[1]; + ray.dir.v[2] = is->vec[2]; + ray.depth = 0; + ray.kdtree_depth = 0; + + /* do this before intersect calls */ + is->vlrcontr= NULL; /* to check shared edge */ + + /* only for shadow! we don't care how far or close; any intersection + * will do (WHY?) */ + if(is->mode==DDA_SHADOW) { + + /* check with last intersected shadow face */ + if(is->vlr_last!=NULL && is->vlr_last!=is->vlrorig) { + if(is->lay & is->vlr_last->lay) { + is->vlr= is->vlr_last; + VECSUB(is->vec, is->end, is->start); + if(intersection(is)) return 1; + } + } + } + + VlakRen* hit; + hit = (VlakRen*) tie_work(&R.tie, &ray, &id, tie_hitfunc, is); + + if (!hit) + /* no intersections found */ + is->vlr_last= NULL; + + return 0; +} + static void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr) { diff -r -u -N -x linux2-config.py -x CVS blender-cvs/source/blender/render/intern/source/tie.c blender/source/blender/render/intern/source/tie.c --- blender-cvs/source/blender/render/intern/source/tie.c 1969-12-31 19:00:00.000000000 -0500 +++ blender/source/blender/render/intern/source/tie.c 2006-06-16 03:38:28.000000000 -0400 @@ -0,0 +1,459 @@ +/* T I E . C +* +* @file tie.c +* +* BRL-CAD +* +* Copyright (c) 2002-2006 United States Government as represented by +* the U.S. Army Research Laboratory. +* +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public License +* as published by the Free Software Foundation; either version 2 of +* the License, or (at your option) any later version. +* +* This library is distributed in the hope that it will be useful, but +* WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +* Library General Public License for more details. +* +* You should have received a copy of the GNU Lesser General Public +* License along with this file; see the file named COPYING for more +* information. +* +* @brief support routines for shooting at triangles +* Comments - +* Triangle Intersection Engine +* +* The calling sequence is as follows: +* - tie_init() initialize the data structure +* - tie_push() add triangles to the universe to be raytraced +* - tie_prep() build the KDTREE for the triangles +* - tie_work() shoot some ray +* - tie_work() shoot some ray +* - tie_work() shoot some ray +* - tie_work() shoot some ray +* - tie_free() Free up all the memory +* +* Author - +* Justin L. Shumaker +* +* Source - +* The U. S. Army Research Laboratory +* Aberdeen Proving Ground, Maryland 21005-5068 USA +* +* $Id: tie.c,v 1.42 2006/01/18 06:46:11 brlcad Exp $ +*/ +/** @addtogroup libtie + * + * The calling sequence is as follows: + * @li tie_init() initialize the data structure + * @li tie_push() add triangles to the universe to be raytraced + * @li tie_prep() build the KDTREE for the triangles + * @li tie_work() shoot some ray + * @li tie_work() shoot some ray + * @li tie_work() shoot some ray + * @li tie_work() shoot some ray + * @li tie_work() shoot some ray + * @li tie_free() Free up all the memory + * @{ + */ + +#ifdef HAVE_CONFIG_H +# include "brlcad_config.h" +#endif + +#include "tie.h" +#include +#include +#include +#include +#include +#include +#include "struct.h" + +#include "MEM_guardedalloc.h" + +#ifdef TIE_SSE +# include +#endif + +/************************************************************* + **************** PRIVATE FUNCTIONS ************************** + *************************************************************/ + +static void tie_tri_prep(tie_t *tie) +{ + TIE_3 v1, v2, u, v; + int i, i1, i2; + tie_tri_t *tri; + + for (i = 0; i < tie->tri_num; i++) { + tri = &tie->tri_list[i]; + + v1 = tri->data[1]; + v2 = tri->data[2]; + + /* Compute Normal */ + math_vec_sub(u, tri->data[1], tri->data[0]); + math_vec_sub(v, tri->data[2], tri->data[0]); + math_vec_cross(tri->data[1], u, v); + math_vec_unitize(tri->data[1]); + + /* Compute i1 and i2 */ + u.v[0] = fabs(tri->data[1].v[0]); + u.v[1] = fabs(tri->data[1].v[1]); + u.v[2] = fabs(tri->data[1].v[2]); + + if (u.v[2] > u.v[1] && u.v[2] > u.v[0]) { + i1 = 0; + i2 = 1; + } else if (u.v[1] > u.v[2] && u.v[1] > u.v[0]) { + i1 = 0; + i2 = 2; + } else { + i1 = 1; + i2 = 2; + } + + /* compute u1, v2, u2, v2 */ + tri->data[2].v[1] = v1.v[i1] - tri->data[0].v[i1]; + tri->data[2].v[2] = v2.v[i1] - tri->data[0].v[i1]; + tri->v[0] = v1.v[i2] - tri->data[0].v[i2]; + tri->v[1] = v2.v[i2] - tri->data[0].v[i2]; + + if (i1 == 0 && i2 == 1) { + tri->v = (tfloat *)((intptr_t)(tri->v) + 2); + } else if (i1 == 0) { + tri->v = (tfloat *)((intptr_t)(tri->v) + 1); + } + + /* Compute DotVN */ + math_vec_mul_scalar(v1, tri->data[0], -1.0); + math_vec_dot(tri->data[2].v[0], v1, tri->data[1]); + } +} + + +/************************************************************* + **************** EXPORTED FUNCTIONS ************************* + *************************************************************/ + +/** + * Initialize a tie_t data structure + * + * This needs to be called before any other libtie data structures are called. + * + * @param tie pointer to a struct tie_t + * @return void + */ +void tie_init(tie_t *tie, unsigned int tri_num) +{ + tie->kdtree = NULL; + tie->tri_num = 0; + tie->tri_num_alloc = tri_num; + tie->tri_list = (tie_tri_t *)MEM_mallocN(sizeof(tie_tri_t) * tri_num, "tri list"); + tie->stat = 0; + tie->rays_fired = 0; +} + +/** + * Free up all the stuff associated with libtie + * + * All of the KDTREE nodes and triangles that we have allocated need to + * be freed in a controlled manner. This routine does that. + * + * @param tie pointer to a struct tie_t + * @return void + */ +void tie_free(tie_t *tie) +{ + int i; + + /* Free Triangle Data */ + for (i = 0; i < tie->tri_num; i++) + MEM_freeN( (void *)((intptr_t)(tie->tri_list[i].v) & ~0x7L)); + MEM_freeN(tie->tri_list); + + /* Free KDTREE Nodes */ + tie_kdtree_free(tie); +} + + +/** + * Get ready to shoot rays at triangles + * + * Build the KDTREE tree for the triangles we have + * + * @param tie pointer to a struct tie_t which now has all the triangles in it + * @return void + */ +void tie_prep(tie_t *tie) +{ + /* Build the kd-tree */ + tie_kdtree_prep(tie); + + /* Prep all the triangles */ + tie_tri_prep(tie); +} + + +/** + * Shoot a ray at some triangles + * + * The user-provided hitfunc is called at each ray/triangle intersection. + * Calls are guaranteed to be made in the ray-intersection order. + * The last argument (void *ptr) is passed to the hitfunc as-is, to allow + * application specific data to be passed to the hitfunc. + * + * @param tie a struct tie_t universe + * @param ray the ray to be intersected with the geometry + * @param id the intersection data for each intersection + * @param hitfunc the application routine to be called upon ray/triangle intersection. + * This function should return 0 if the ray is to continue propagating through the geometry, + * or non-zero if ray intersection should cease. + * @param ptr a pointer to be passed to the hitfunc when it is called. + * + * @return the return value from the user hitfunc() is used. + * In the event that the ray does not hit anything, or the ray exits the geometry space, a null value will be returned. + * @retval 0 ray did not hit anything, or ray was propagated through the geometry completely. + * @retval !0 the value returned from the last invokation of hitfunc() + */ +void* tie_work(tie_t *tie, tie_ray_t *ray, tie_id_t *id, void *(*hitfunc)(tie_ray_t*, tie_id_t*, tie_tri_t*, void *ptr), void *ptr) +{ + tie_stack_t stack[40]; + tie_id_t t, id_list[256]; + tie_tri_t *hit_list[256], *tri; + tie_geom_t *data; + tie_kdtree_t *node_aligned, *temp[2]; + tfloat near, far, dirinv[3], dist; + int i, n, ab[3], split, stack_ind, hit_count; + void *result; + + + if (!tie->kdtree) + return (NULL); + + ray->kdtree_depth = 0; + + /* + * Precompute direction inverse since it's used in a bunch of divides, + * this allows those divides to become fast multiplies. + */ + for (i = 0; i < 3; i++) { + if (ray->dir.v[i] == 0) + ray->dir.v[i] = TIE_PREC; + dirinv[i] = 1.0 / ray->dir.v[i]; + ab[i] = dirinv[i] < 0 ? 1 : 0; + } + + /* Extracting value of splitting plane from tie->kdtree pointer */ + split = ((intptr_t)(((tie_kdtree_t *)((intptr_t)tie->kdtree & ~0x7L))->data)) & 0x3; + + /* Initialize ray segment */ + if (ray->dir.v[split] < 0) { + far = (tie->min.v[split] - ray->pos.v[split]) * dirinv[split]; + } else { + far = (tie->max.v[split] - ray->pos.v[split]) * dirinv[split]; + } + + stack_ind = 0; + stack[0].node = tie->kdtree; + stack[0].near = 0; + stack[0].far = far; + + /* Process items on the stack */ + while (stack_ind >= 0) { + near = stack[stack_ind].near; + far = stack[stack_ind].far; + + /* + * Take the pointer from stack[stack_ind] and remove lower pts bits used to store data to + * give a valid ptr address. + */ + node_aligned = (tie_kdtree_t *)((intptr_t)stack[stack_ind].node & ~0x7L); + stack_ind--; + + /* + * KDTREE TRAVERSAL + * + * 3 conditions can happen here: + * - Ray only intersects the nearest node + * - Ray only intersects the furthest node + * - Ray intersects both nodes, pushing the furthest onto the stack + * + * Gordon Stoll's Mantra - Rays are Measured in Millions :-) + */ + while (((intptr_t)(node_aligned->data)) & 0x4) { + ray->kdtree_depth++; + + /* Retreive the splitting plane */ + split = ((intptr_t)(node_aligned->data)) & 0x3; + + /* Calculate the projected 1d distance to splitting axis */ + dist = (node_aligned->axis - ray->pos.v[split]) * dirinv[split]; + + temp[0] = &((tie_kdtree_t *)(node_aligned->data))[ab[split]]; + temp[1] = &((tie_kdtree_t *)(node_aligned->data))[1 - ab[split]]; + + i = near >= dist; /* Node B Only? */ + node_aligned = (tie_kdtree_t *)((intptr_t)(temp[i]) & ~0x7L); + + if (far < dist || i) + continue; + + /* Nearest Node and Push Furthest */ + stack_ind++; + stack[stack_ind].node = temp[1]; + stack[stack_ind].near = dist; + stack[stack_ind].far = far; + far = dist; + } + + + /* + * RAY/TRIANGLE INTERSECTION - Only gets executed on geometry nodes. + * This part of the function is being executed because the KDTREE Traversal is Complete. + */ + data = (tie_geom_t *)(node_aligned->data); + if (data->tri_num == 0) + continue; + + hit_count = 0; + for (i = 0; i < data->tri_num; i++) { + /* + * Triangle Intersection Code + */ + tfloat u0, v0, *v; + int i1, i2; + + tri = data->tri_list[i]; + math_vec_dot(u0, tri->data[1], ray->pos); + math_vec_dot(v0, tri->data[1], ray->dir); + t.dist = -(tri->data[2].v[0] + u0) / v0; + + /* + * Intersection point on triangle must lie within the kdtree node or it is rejected + * Apply TIE_PREC to near and far such that triangles that lie on orthogonal planes + * aren't in a precision fuzz boundary, thus missing something they should actualy + * have hit. + */ + if (t.dist < near - TIE_PREC || t.dist > far + TIE_PREC) + continue; + + /* Compute Intersection Point (P = O + Dt) */ + math_vec_mul_scalar(t.pos, ray->dir, t.dist); + math_vec_add(t.pos, ray->pos, t.pos); + + /* Extract i1 and i2 indices from lower bits of the v pointer */ + v = (tfloat *)((intptr_t)(tri->v) & ~0x7L); + + i1 = TIE_TAB1[((intptr_t)(tri->v) & 0x7)]; + i2 = TIE_TAB1[3 + ((intptr_t)(tri->v) & 0x7)]; + + /* Compute U and V */ + u0 = t.pos.v[i1] - tri->data[0].v[i1]; + v0 = t.pos.v[i2] - tri->data[0].v[i2]; + + /* + * Compute the barycentric coordinates, and make sure the coordinates + * fall within the boundaries of the triangle plane. + */ + if (fabs(tri->data[2].v[1]) <= TIE_PREC) { + t.beta = u0 / tri->data[2].v[2]; + if (t.beta < 0 || t.beta > 1) + continue; + t.alpha = (v0 - t.beta * v[1]) / v[0]; + } else { + t.beta = (v0 * tri->data[2].v[1] - u0 * v[0]) / (v[1] * tri->data[2].v[1] - tri->data[2].v[2] * v[0]); + if (t.beta < 0 || t.beta > 1) + continue; + t.alpha = (u0 - t.beta * tri->data[2].v[2]) / tri->data[2].v[1]; + } + + if (t.alpha < 0 || (t.alpha + t.beta) > 1.0) + continue; + + /* Triangle Intersected, append it in the list */ + if (hit_count < 0xff) { + hit_list[hit_count] = tri; + id_list[hit_count] = t; + hit_count++; + } + + + } + + + /* If we hit something, then sort the hit triangles on demand */ + for (i = 0; i < hit_count; i++) { + /* Sort the list so that HitList and IDList [n] is in order wrt [i] */ + for (n = i; n < hit_count; n++) { + if (id_list[n].dist < id_list[i].dist) { + /* Swap */ + tri = hit_list[i]; + t = id_list[i]; + hit_list[i] = hit_list[n]; + id_list[i] = id_list[n]; + hit_list[n] = tri; + id_list[n] = t; + } + } + + id_list[i].norm = hit_list[i]->data[1]; + result = hitfunc(ray, &id_list[i], hit_list[i], ptr); + + if (result) { + *id = id_list[i]; + return (result); + } + } + } + + return (NULL); +} + + +/** + * Add a triangle + * + * Add a new triangle to the universe to be raytraced. + * + * @param tie the universe + * @param tlist is an array of TIE_3 vertice triplets (v0, v1, v2) that form each triangle. + * @param tnum is the number of triangles (tlist = 3 * tnum of TIE_3's). + * @param plist is a list of pointer data that gets assigned to the ptr of each triangle. + * This will typically be 4-byte (32-bit architecture) spaced array of pointers that + * associate the triangle pointer with your arbitrary structure, i.e a mesh. + * @param pstride is the number of bytes to increment the pointer list as it assigns + * a pointer to each mesh, typically a value of 4 (for 32-bit machines). If you have + * a single pointer that groups all triangles to a common structure then you can use + * a value of 0 for pstride. This will give the pointer of all triangles the pointer + * address of plist. + * @return void + */ +void tie_push(tie_t *tie, TIE_3 *tlist, int tnum, void *plist, int pstride) +{ + int i; + + if (tnum + tie->tri_num > tie->tri_num_alloc) { + tie->tri_list = (tie_tri_t *)MEM_reallocN(tie->tri_list, sizeof(tie_tri_t) * (tie->tri_num + tnum)); + tie->tri_num_alloc += tnum; + } + + for (i = 0; i < tnum; i++) { + tie->tri_list[tie->tri_num].data[0] = tlist[i * 3 + 0]; + tie->tri_list[tie->tri_num].data[1] = tlist[i * 3 + 1]; + tie->tri_list[tie->tri_num].data[2] = tlist[i * 3 + 2]; + if (plist) { + tie->tri_list[tie->tri_num].ptr = plist; + plist = (void *)((intptr_t)plist + pstride); + } else { + tie->tri_list[tie->tri_num].ptr = NULL; + } + tie->tri_list[tie->tri_num].v = (tfloat *)MEM_mallocN(2 * sizeof(tfloat), "kdtree tri_list push"); + tie->tri_num++; + } +} + +/** @} */