clray
changeset 54:6a30f27fa1e6
separated the OpenGL visualization and added a CPU raytracing mode
| author | John Tsiombikas <nuclear@member.fsf.org> |
|---|---|
| date | Fri, 10 Sep 2010 16:47:00 +0100 |
| parents | 54a96b738afe |
| children | df239a52a091 |
| files | rt.cl src/clray.cc src/common.h src/dbggl.cc src/dbgray.cc src/rt.cc src/rt.cl src/rt.h src/scene.cc src/scene.h src/vector.cc src/vector.h src/vector.inl |
| diffstat | 13 files changed, 977 insertions(+), 498 deletions(-) [+] |
line diff
1.1 --- a/rt.cl Sun Sep 05 16:43:55 2010 +0100 1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 @@ -1,395 +0,0 @@ 1.4 -/* vim: set ft=opencl:ts=4:sw=4 */ 1.5 -#include "common.h" 1.6 - 1.7 -struct RendInfo { 1.8 - float4 ambient; 1.9 - int xsz, ysz; 1.10 - int num_faces, num_lights; 1.11 - int max_iter; 1.12 - int cast_shadows; 1.13 -}; 1.14 - 1.15 -struct Vertex { 1.16 - float4 pos; 1.17 - float4 normal; 1.18 - float4 tex; 1.19 - float4 padding; 1.20 -}; 1.21 - 1.22 -struct Face { 1.23 - struct Vertex v[3]; 1.24 - float4 normal; 1.25 - int matid; 1.26 - int padding[3]; 1.27 -}; 1.28 - 1.29 -struct Material { 1.30 - float4 kd, ks; 1.31 - float kr, kt; 1.32 - float spow; 1.33 - float padding; 1.34 -}; 1.35 - 1.36 -struct Light { 1.37 - float4 pos, color; 1.38 -}; 1.39 - 1.40 -struct Ray { 1.41 - float4 origin, dir; 1.42 -}; 1.43 - 1.44 -struct SurfPoint { 1.45 - float t; 1.46 - float4 pos, norm, dbg; 1.47 - global const struct Face *obj; 1.48 - struct Material mat; 1.49 -}; 1.50 - 1.51 -struct Scene { 1.52 - float4 ambient; 1.53 - global const struct Face *faces; 1.54 - int num_faces; 1.55 - global const struct Light *lights; 1.56 - int num_lights; 1.57 - global const struct Material *matlib; 1.58 - //global const struct KDNode *kdtree; 1.59 - bool cast_shadows; 1.60 -}; 1.61 - 1.62 -struct AABBox { 1.63 - float4 min, max; 1.64 -}; 1.65 - 1.66 -struct KDNode { 1.67 - struct AABBox aabb; 1.68 - int face_idx[MAX_NODE_FACES]; 1.69 - int num_faces; 1.70 - int left, right; 1.71 - int padding; 1.72 -}; 1.73 - 1.74 -#define MIN_ENERGY 0.001 1.75 -#define EPSILON 1e-5 1.76 - 1.77 -float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg); 1.78 -bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *sp, read_only image2d_t kdimg); 1.79 -bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp); 1.80 -bool intersect_aabb(struct Ray ray, struct AABBox aabb); 1.81 - 1.82 -float4 reflect(float4 v, float4 n); 1.83 -float4 transform(float4 v, global const float *xform); 1.84 -void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans); 1.85 -float4 calc_bary(float4 pt, global const struct Face *face, float4 norm); 1.86 -float mean(float4 v); 1.87 - 1.88 -void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg); 1.89 - 1.90 - 1.91 -kernel void render(write_only image2d_t fb, 1.92 - global const struct RendInfo *rinf, 1.93 - global const struct Face *faces, 1.94 - global const struct Material *matlib, 1.95 - global const struct Light *lights, 1.96 - global const struct Ray *primrays, 1.97 - global const float *xform, 1.98 - global const float *invtrans, 1.99 - //global const struct KDNode *kdtree 1.100 - read_only image2d_t kdtree_img) 1.101 -{ 1.102 - int idx = get_global_id(0); 1.103 - 1.104 - struct Scene scn; 1.105 - scn.ambient = rinf->ambient; 1.106 - scn.faces = faces; 1.107 - scn.num_faces = rinf->num_faces; 1.108 - scn.lights = lights; 1.109 - scn.num_lights = rinf->num_lights; 1.110 - scn.matlib = matlib; 1.111 - scn.cast_shadows = rinf->cast_shadows; 1.112 - 1.113 - struct Ray ray = primrays[idx]; 1.114 - transform_ray(&ray, xform, invtrans); 1.115 - 1.116 - float4 pixel = (float4)(0, 0, 0, 0); 1.117 - float4 energy = (float4)(1.0, 1.0, 1.0, 0.0); 1.118 - int iter = 0; 1.119 - 1.120 - while(iter++ <= rinf->max_iter && mean(energy) > MIN_ENERGY) { 1.121 - struct SurfPoint sp; 1.122 - if(find_intersection(ray, &scn, &sp, kdtree_img)) { 1.123 - pixel += shade(ray, &scn, &sp, kdtree_img) * energy; 1.124 - 1.125 - float4 refl_col = sp.mat.ks * sp.mat.kr; 1.126 - 1.127 - ray.origin = sp.pos; 1.128 - ray.dir = reflect(-ray.dir, sp.norm); 1.129 - 1.130 - energy *= refl_col; 1.131 - } else { 1.132 - energy = (float4)(0.0, 0.0, 0.0, 0.0); 1.133 - } 1.134 - } 1.135 - 1.136 - int2 coord; 1.137 - coord.x = idx % rinf->xsz; 1.138 - coord.y = idx / rinf->xsz; 1.139 - 1.140 - write_imagef(fb, coord, pixel); 1.141 -} 1.142 - 1.143 -float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg) 1.144 -{ 1.145 - float4 norm = sp->norm; 1.146 - //bool entering = true; 1.147 - 1.148 - if(dot(ray.dir, norm) >= 0.0) { 1.149 - norm = -norm; 1.150 - //entering = false; 1.151 - } 1.152 - 1.153 - float4 dcol = scn->ambient * sp->mat.kd; 1.154 - float4 scol = (float4)(0, 0, 0, 0); 1.155 - 1.156 - for(int i=0; i<scn->num_lights; i++) { 1.157 - float4 ldir = scn->lights[i].pos - sp->pos; 1.158 - 1.159 - struct Ray shadowray; 1.160 - shadowray.origin = sp->pos; 1.161 - shadowray.dir = ldir; 1.162 - 1.163 - if(!scn->cast_shadows || !find_intersection(shadowray, scn, 0, kdimg)) { 1.164 - ldir = normalize(ldir); 1.165 - float4 vdir = -ray.dir; 1.166 - vdir.x = native_divide(vdir.x, RAY_MAG); 1.167 - vdir.y = native_divide(vdir.y, RAY_MAG); 1.168 - vdir.z = native_divide(vdir.z, RAY_MAG); 1.169 - float4 vref = reflect(vdir, norm); 1.170 - 1.171 - float diff = fmax(dot(ldir, norm), 0.0f); 1.172 - dcol += sp->mat.kd /* scn->lights[i].color*/ * diff; 1.173 - 1.174 - float spec = native_powr(fmax(dot(ldir, vref), 0.0f), sp->mat.spow); 1.175 - scol += sp->mat.ks /* scn->lights[i].color*/ * spec; 1.176 - } 1.177 - } 1.178 - 1.179 - return dcol + scol; 1.180 -} 1.181 - 1.182 -#define STACK_SIZE MAX_TREE_DEPTH 1.183 -bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *spres, read_only image2d_t kdimg) 1.184 -{ 1.185 - struct SurfPoint sp0; 1.186 - sp0.t = 1.0; 1.187 - sp0.obj = 0; 1.188 - 1.189 - int idxstack[STACK_SIZE]; 1.190 - int top = 0; // points after the topmost element of the stack 1.191 - idxstack[top++] = 0; // root at tree[0] 1.192 - 1.193 - while(top > 0) { 1.194 - int idx = idxstack[--top]; // remove this index from the stack and process it 1.195 - 1.196 - struct KDNode node; 1.197 - read_kdnode(idx, &node, kdimg); 1.198 - 1.199 - if(intersect_aabb(ray, node.aabb)) { 1.200 - if(node.left == -1) { 1.201 - // leaf node... check each face in turn and update the nearest intersection as needed 1.202 - for(int i=0; i<node.num_faces; i++) { 1.203 - struct SurfPoint spt; 1.204 - int fidx = node.face_idx[i]; 1.205 - 1.206 - if(intersect(ray, scn->faces + fidx, &spt) && spt.t < sp0.t) { 1.207 - sp0 = spt; 1.208 - } 1.209 - } 1.210 - } else { 1.211 - // internal node... recurse to the children 1.212 - idxstack[top++] = node.left; 1.213 - idxstack[top++] = node.right; 1.214 - } 1.215 - } 1.216 - } 1.217 - 1.218 - if(!sp0.obj) { 1.219 - return false; 1.220 - } 1.221 - 1.222 - if(spres) { 1.223 - *spres = sp0; 1.224 - spres->mat = scn->matlib[sp0.obj->matid]; 1.225 - } 1.226 - return true; 1.227 -} 1.228 - 1.229 -bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp) 1.230 -{ 1.231 - float4 origin = ray.origin; 1.232 - float4 dir = ray.dir; 1.233 - float4 norm = face->normal; 1.234 - 1.235 - float ndotdir = dot(dir, norm); 1.236 - 1.237 - if(fabs(ndotdir) <= EPSILON) { 1.238 - return false; 1.239 - } 1.240 - 1.241 - float4 pt = face->v[0].pos; 1.242 - float4 vec = pt - origin; 1.243 - 1.244 - float ndotvec = dot(norm, vec); 1.245 - float t = native_divide(ndotvec, ndotdir); 1.246 - 1.247 - if(t < EPSILON || t > 1.0) { 1.248 - return false; 1.249 - } 1.250 - pt = origin + dir * t; 1.251 - 1.252 - 1.253 - float4 bc = calc_bary(pt, face, norm); 1.254 - float bc_sum = bc.x + bc.y + bc.z; 1.255 - 1.256 - if(bc_sum < 1.0 - EPSILON || bc_sum > 1.0 + EPSILON) { 1.257 - return false; 1.258 - bc *= 1.2; 1.259 - } 1.260 - 1.261 - sp->t = t; 1.262 - sp->pos = pt; 1.263 - sp->norm = normalize(face->v[0].normal * bc.x + face->v[1].normal * bc.y + face->v[2].normal * bc.z); 1.264 - sp->obj = face; 1.265 - sp->dbg = bc; 1.266 - return true; 1.267 -} 1.268 - 1.269 -bool intersect_aabb(struct Ray ray, struct AABBox aabb) 1.270 -{ 1.271 - if(ray.origin.x >= aabb.min.x && ray.origin.y >= aabb.min.y && ray.origin.z >= aabb.min.z && 1.272 - ray.origin.x < aabb.max.x && ray.origin.y < aabb.max.y && ray.origin.z < aabb.max.z) { 1.273 - return true; 1.274 - } 1.275 - 1.276 - float4 bbox[2] = { 1.277 - aabb.min.x, aabb.min.y, aabb.min.z, 0, 1.278 - aabb.max.x, aabb.max.y, aabb.max.z, 0 1.279 - }; 1.280 - 1.281 - int xsign = (int)(ray.dir.x < 0.0); 1.282 - float invdirx = native_recip(ray.dir.x); 1.283 - float tmin = (bbox[xsign].x - ray.origin.x) * invdirx; 1.284 - float tmax = (bbox[1 - xsign].x - ray.origin.x) * invdirx; 1.285 - 1.286 - int ysign = (int)(ray.dir.y < 0.0); 1.287 - float invdiry = native_recip(ray.dir.y); 1.288 - float tymin = (bbox[ysign].y - ray.origin.y) * invdiry; 1.289 - float tymax = (bbox[1 - ysign].y - ray.origin.y) * invdiry; 1.290 - 1.291 - if(tmin > tymax || tymin > tmax) { 1.292 - return false; 1.293 - } 1.294 - 1.295 - if(tymin > tmin) tmin = tymin; 1.296 - if(tymax < tmax) tmax = tymax; 1.297 - 1.298 - int zsign = (int)(ray.dir.z < 0.0); 1.299 - float invdirz = native_recip(ray.dir.z); 1.300 - float tzmin = (bbox[zsign].z - ray.origin.z) * invdirz; 1.301 - float tzmax = (bbox[1 - zsign].z - ray.origin.z) * invdirz; 1.302 - 1.303 - if(tmin > tzmax || tzmin > tmax) { 1.304 - return false; 1.305 - } 1.306 - 1.307 - return tmin < 1.0 && tmax > 0.0; 1.308 -} 1.309 - 1.310 -float4 reflect(float4 v, float4 n) 1.311 -{ 1.312 - return 2.0f * dot(v, n) * n - v; 1.313 -} 1.314 - 1.315 -float4 transform(float4 v, global const float *xform) 1.316 -{ 1.317 - float4 res; 1.318 - res.x = v.x * xform[0] + v.y * xform[4] + v.z * xform[8] + xform[12]; 1.319 - res.y = v.x * xform[1] + v.y * xform[5] + v.z * xform[9] + xform[13]; 1.320 - res.z = v.x * xform[2] + v.y * xform[6] + v.z * xform[10] + xform[14]; 1.321 - res.w = 0.0; 1.322 - return res; 1.323 -} 1.324 - 1.325 -void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans) 1.326 -{ 1.327 - ray->origin = transform(ray->origin, xform); 1.328 - ray->dir = transform(ray->dir, invtrans); 1.329 -} 1.330 - 1.331 -float4 calc_bary(float4 pt, global const struct Face *face, float4 norm) 1.332 -{ 1.333 - float4 bc = (float4)(0, 0, 0, 0); 1.334 - 1.335 - // calculate area of the whole triangle 1.336 - float4 v1 = face->v[1].pos - face->v[0].pos; 1.337 - float4 v2 = face->v[2].pos - face->v[0].pos; 1.338 - float4 xv1v2 = cross(v1, v2); 1.339 - 1.340 - float area = fabs(dot(xv1v2, norm)) * 0.5; 1.341 - if(area < EPSILON) { 1.342 - return bc; 1.343 - } 1.344 - 1.345 - float4 pv0 = face->v[0].pos - pt; 1.346 - float4 pv1 = face->v[1].pos - pt; 1.347 - float4 pv2 = face->v[2].pos - pt; 1.348 - 1.349 - // calculate the area of each sub-triangle 1.350 - float4 x12 = cross(pv1, pv2); 1.351 - float4 x20 = cross(pv2, pv0); 1.352 - float4 x01 = cross(pv0, pv1); 1.353 - 1.354 - float a0 = fabs(dot(x12, norm)) * 0.5; 1.355 - float a1 = fabs(dot(x20, norm)) * 0.5; 1.356 - float a2 = fabs(dot(x01, norm)) * 0.5; 1.357 - 1.358 - bc.x = native_divide(a0, area); 1.359 - bc.y = native_divide(a1, area); 1.360 - bc.z = native_divide(a2, area); 1.361 - return bc; 1.362 -} 1.363 - 1.364 -float mean(float4 v) 1.365 -{ 1.366 - return native_divide(v.x + v.y + v.z, 3.0); 1.367 -} 1.368 - 1.369 - 1.370 -const sampler_t kdsampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_NONE | CLK_FILTER_NEAREST; 1.371 - 1.372 -// read a KD-tree node from a texture scanline 1.373 -void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg) 1.374 -{ 1.375 - int startx = KDIMG_NODE_WIDTH * (idx / KDIMG_MAX_HEIGHT); 1.376 - 1.377 - int2 tc; 1.378 - tc.x = startx; 1.379 - tc.y = idx % KDIMG_MAX_HEIGHT; 1.380 - 1.381 - node->aabb.min = read_imagef(kdimg, kdsampler, tc); tc.x++; 1.382 - node->aabb.max = read_imagef(kdimg, kdsampler, tc); 1.383 - 1.384 - tc.x = startx + 2 + MAX_NODE_FACES / 4; 1.385 - float4 pix = read_imagef(kdimg, kdsampler, tc); 1.386 - node->num_faces = (int)pix.x; 1.387 - node->left = (int)pix.y; 1.388 - node->right = (int)pix.z; 1.389 - 1.390 - tc.x = startx + 2; 1.391 - for(int i=0; i<node->num_faces; i+=4) { 1.392 - float4 pix = read_imagef(kdimg, kdsampler, tc); tc.x++; 1.393 - node->face_idx[i] = (int)pix.x; 1.394 - node->face_idx[i + 1] = (int)pix.y; 1.395 - node->face_idx[i + 2] = (int)pix.z; 1.396 - node->face_idx[i + 3] = (int)pix.w; 1.397 - } 1.398 -}
2.1 --- a/src/clray.cc Sun Sep 05 16:43:55 2010 +0100 2.2 +++ b/src/clray.cc Fri Sep 10 16:47:00 2010 +0100 2.3 @@ -33,14 +33,21 @@ 2.4 static float cam_theta, cam_phi = 25.0; 2.5 static float cam_dist = 10.0; 2.6 2.7 -static bool dbg_glrender = false; 2.8 -static bool dbg_show_kdtree = false; 2.9 +static bool dbg_glrender; 2.10 +static bool dbg_nocl; 2.11 +static bool dbg_show_kdtree; 2.12 static bool dbg_show_obj = true; 2.13 bool dbg_frame_time = true; 2.14 2.15 static Scene scn; 2.16 static unsigned int tex; 2.17 2.18 +static Light lightlist[] = { 2.19 + {{-8, 15, 18, 0}, {1, 1, 1, 1}} 2.20 +}; 2.21 + 2.22 + 2.23 + 2.24 int main(int argc, char **argv) 2.25 { 2.26 glutInitWindowSize(800, 600); 2.27 @@ -81,6 +88,10 @@ 2.28 dbg_glrender = true; 2.29 break; 2.30 2.31 + case 'n': 2.32 + dbg_nocl = true; 2.33 + break; 2.34 + 2.35 default: 2.36 fprintf(stderr, "unrecognized option: %s\n", argv[i]); 2.37 return 1; 2.38 @@ -103,6 +114,11 @@ 2.39 return false; 2.40 } 2.41 2.42 + int num_lights = sizeof lightlist / sizeof *lightlist; 2.43 + for(int i=0; i<num_lights; i++) { 2.44 + scn.add_light(lightlist[i]); 2.45 + } 2.46 + 2.47 glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE); 2.48 glutCreateWindow("OpenCL Raytracer"); 2.49 2.50 @@ -171,11 +187,6 @@ 2.51 2.52 glGetFloatv(GL_MODELVIEW_MATRIX, mat.m); 2.53 2.54 - inv_mat = mat; 2.55 - inv_mat.invert(); 2.56 - 2.57 - /*inv_trans = inv_mat; 2.58 - inv_trans.transpose();*/ 2.59 inv_trans = mat; 2.60 inv_trans.m[3] = inv_trans.m[7] = inv_trans.m[11] = 0.0; 2.61 inv_trans.m[12] = inv_trans.m[13] = inv_trans.m[14] = 0.0; 2.62 @@ -185,14 +196,21 @@ 2.63 glPopMatrix(); 2.64 2.65 if(!dbg_glrender) { 2.66 - if(!render()) { 2.67 - exit(1); 2.68 + if(dbg_nocl) { 2.69 + dbg_render(mat.m, inv_trans.m); 2.70 + } else { 2.71 + if(!render()) { 2.72 + exit(1); 2.73 + } 2.74 } 2.75 need_update = false; 2.76 } 2.77 } 2.78 2.79 if(dbg_glrender) { 2.80 + inv_mat = mat; 2.81 + inv_mat.invert(); 2.82 + 2.83 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); 2.84 glLoadMatrixf(inv_mat.m); 2.85 dbg_render_gl(&scn, dbg_show_kdtree, dbg_show_obj); 2.86 @@ -322,6 +340,13 @@ 2.87 dbg_frame_time = !dbg_frame_time; 2.88 break; 2.89 2.90 + case 'n': 2.91 + dbg_nocl = !dbg_nocl; 2.92 + printf("switching to %s rendering\n", dbg_nocl ? "debug CPU" : "OpenCL"); 2.93 + need_update = true; 2.94 + glutPostRedisplay(); 2.95 + break; 2.96 + 2.97 default: 2.98 break; 2.99 }
3.1 --- a/src/common.h Sun Sep 05 16:43:55 2010 +0100 3.2 +++ b/src/common.h Fri Sep 10 16:47:00 2010 +0100 3.3 @@ -1,6 +1,12 @@ 3.4 #ifndef COMMON_H_ 3.5 #define COMMON_H_ 3.6 3.7 +/* error threshold */ 3.8 +#define EPSILON 1e-5 3.9 + 3.10 +/* minimum trace energy threshold */ 3.11 +#define MIN_ENERGY 0.001 3.12 + 3.13 /* primary ray magnitude */ 3.14 #define RAY_MAG 500.0 3.15
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 4.2 +++ b/src/dbggl.cc Fri Sep 10 16:47:00 2010 +0100 4.3 @@ -0,0 +1,77 @@ 4.4 +#include <string.h> 4.5 +#include <assert.h> 4.6 +#include "rt.h" 4.7 +#include "ogl.h" 4.8 + 4.9 +#define MIN(a, b) ((a) < (b) ? (a) : (b)) 4.10 +static void dbg_set_gl_material(Material *mat) 4.11 +{ 4.12 + static Material def_mat = {{0.7, 0.7, 0.7, 1}, {0, 0, 0, 0}, 0, 0, 0}; 4.13 + 4.14 + if(!mat) mat = &def_mat; 4.15 + 4.16 + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat->kd); 4.17 + glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat->ks); 4.18 + glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, MIN(mat->spow, 128.0f)); 4.19 +} 4.20 + 4.21 +void dbg_render_gl(Scene *scn, bool show_tree, bool show_obj) 4.22 +{ 4.23 + const RendInfo *rinf = get_render_info(); 4.24 + const Face *faces = scn->get_face_buffer(); 4.25 + 4.26 + glPushAttrib(GL_ENABLE_BIT | GL_TRANSFORM_BIT | GL_LIGHTING_BIT); 4.27 + 4.28 + for(int i=0; i<scn->get_num_lights(); i++) { 4.29 + float lpos[4]; 4.30 + 4.31 + memcpy(lpos, scn->lights[i].pos, sizeof lpos); 4.32 + lpos[3] = 1.0; 4.33 + 4.34 + glLightfv(GL_LIGHT0 + i, GL_POSITION, lpos); 4.35 + glLightfv(GL_LIGHT0 + i, GL_DIFFUSE, scn->lights[i].color); 4.36 + glEnable(GL_LIGHT0 + i); 4.37 + } 4.38 + 4.39 + glDisable(GL_TEXTURE_2D); 4.40 + glEnable(GL_DEPTH_TEST); 4.41 + glEnable(GL_LIGHTING); 4.42 + 4.43 + glMatrixMode(GL_PROJECTION); 4.44 + glPushMatrix(); 4.45 + glLoadIdentity(); 4.46 + gluPerspective(45.0, (float)rinf->xsz / (float)rinf->ysz, 0.5, 1000.0); 4.47 + 4.48 + if(show_obj) { 4.49 + Material *materials = scn->get_materials(); 4.50 + 4.51 + int num_faces = scn->get_num_faces(); 4.52 + int cur_mat = -1; 4.53 + 4.54 + for(int i=0; i<num_faces; i++) { 4.55 + if(faces[i].matid != cur_mat) { 4.56 + if(cur_mat != -1) { 4.57 + glEnd(); 4.58 + } 4.59 + dbg_set_gl_material(materials ? materials + faces[i].matid : 0); 4.60 + cur_mat = faces[i].matid; 4.61 + glBegin(GL_TRIANGLES); 4.62 + } 4.63 + 4.64 + for(int j=0; j<3; j++) { 4.65 + glNormal3fv(faces[i].v[j].normal); 4.66 + glVertex3fv(faces[i].v[j].pos); 4.67 + } 4.68 + } 4.69 + glEnd(); 4.70 + } 4.71 + 4.72 + if(show_tree) { 4.73 + scn->draw_kdtree(); 4.74 + } 4.75 + 4.76 + glPopMatrix(); 4.77 + glPopAttrib(); 4.78 + 4.79 + assert(glGetError() == GL_NO_ERROR); 4.80 +}
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 5.2 +++ b/src/dbgray.cc Fri Sep 10 16:47:00 2010 +0100 5.3 @@ -0,0 +1,342 @@ 5.4 +#include <string.h> 5.5 +#include <assert.h> 5.6 +#include "rt.h" 5.7 +#include "ogl.h" 5.8 +#include "vector.h" 5.9 +#include "timer.h" 5.10 + 5.11 +struct SurfPoint { 5.12 + float t; 5.13 + Vector3 pos, norm; 5.14 + const Face *face; 5.15 +}; 5.16 + 5.17 +static void trace_ray(float *pixel, const Ray &ray, int iter, float energy = 1.0f); 5.18 +static void shade(float *pixel, const Ray &ray, const SurfPoint &sp, int iter, float energy = 1.0f); 5.19 +static bool find_intersection(const Ray &ray, const Scene *scn, const KDNode *kd, SurfPoint *spret); 5.20 +static bool ray_aabb_test(const Ray &ray, const AABBox &aabb); 5.21 +static bool ray_triangle_test(const Ray &ray, const Face *face, SurfPoint *sp); 5.22 +static Vector3 calc_bary(const Vector3 &pt, const Face *face, const Vector3 &norm); 5.23 +static void transform(float *res, const float *v, const float *xform); 5.24 +static void transform_ray(Ray *ray, const float *xform, const float *invtrans_xform); 5.25 + 5.26 +static int xsz, ysz; 5.27 +static float *fb; 5.28 +static unsigned int tex; 5.29 +static Scene *scn; 5.30 +static const Ray *prim_rays; 5.31 + 5.32 + 5.33 +bool init_dbg_renderer(int width, int height, Scene *scene, unsigned int texid) 5.34 +{ 5.35 + try { 5.36 + fb = new float[3 * width * height]; 5.37 + } 5.38 + catch(...) { 5.39 + return false; 5.40 + } 5.41 + 5.42 + xsz = width; 5.43 + ysz = height; 5.44 + tex = texid; 5.45 + scn = scene; 5.46 + return true; 5.47 +} 5.48 + 5.49 +void destroy_dbg_renderer() 5.50 +{ 5.51 + delete [] fb; 5.52 + delete [] prim_rays; 5.53 +} 5.54 + 5.55 +void dbg_set_primary_rays(const Ray *rays) 5.56 +{ 5.57 + prim_rays = rays; 5.58 +} 5.59 + 5.60 +void dbg_render(const float *xform, const float *invtrans_xform, int num_threads) 5.61 +{ 5.62 + unsigned long t0 = get_msec(); 5.63 + 5.64 + int iter = get_render_option_int(ROPT_ITER); 5.65 + 5.66 + int offs = 0; 5.67 + for(int i=0; i<ysz; i++) { 5.68 + for(int j=0; j<xsz; j++) { 5.69 + Ray ray = prim_rays[offs]; 5.70 + transform_ray(&ray, xform, invtrans_xform); 5.71 + 5.72 + trace_ray(fb + offs * 3, ray, iter, 1.0); 5.73 + offs++; 5.74 + } 5.75 + } 5.76 + 5.77 + glPushAttrib(GL_TEXTURE_BIT); 5.78 + glBindTexture(GL_TEXTURE_2D, tex); 5.79 + glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, xsz, ysz, GL_RGB, GL_FLOAT, fb); 5.80 + glPopAttrib(); 5.81 + 5.82 + printf("rendered in %lu msec\n", get_msec() - t0); 5.83 +} 5.84 + 5.85 +static void trace_ray(float *pixel, const Ray &ray, int iter, float energy) 5.86 +{ 5.87 + SurfPoint sp; 5.88 + 5.89 + if(find_intersection(ray, scn, scn->kdtree, &sp)) { 5.90 + shade(pixel, ray, sp, iter, energy); 5.91 + } else { 5.92 + pixel[0] = pixel[1] = pixel[2] = 0.05f; 5.93 + } 5.94 +} 5.95 + 5.96 +#define MAX(a, b) ((a) > (b) ? (a) : (b)) 5.97 + 5.98 +static void shade(float *pixel, const Ray &ray, const SurfPoint &sp, int iter, float energy) 5.99 +{ 5.100 + const Material *mat = scn->get_materials() + sp.face->matid; 5.101 + 5.102 + bool cast_shadows = get_render_option_bool(ROPT_SHAD); 5.103 + Vector3 raydir(ray.dir); 5.104 + Vector3 norm = sp.norm; 5.105 + 5.106 + if(dot(raydir, norm) >= 0.0) { 5.107 + norm = -norm; 5.108 + } 5.109 + 5.110 + float dcol[3] = {0, 0, 0}; 5.111 + float scol[3] = {0, 0, 0}; 5.112 + 5.113 + for(int i=0; i<scn->get_num_lights(); i++) { 5.114 + Vector3 lpos(scn->lights[i].pos); 5.115 + Vector3 ldir = lpos - sp.pos; 5.116 + 5.117 + Ray shadowray; 5.118 + shadowray.origin[0] = sp.pos.x; 5.119 + shadowray.origin[1] = sp.pos.y; 5.120 + shadowray.origin[2] = sp.pos.z; 5.121 + shadowray.dir[0] = ldir.x; 5.122 + shadowray.dir[1] = ldir.y; 5.123 + shadowray.dir[2] = ldir.z; 5.124 + 5.125 + if(!cast_shadows || !find_intersection(shadowray, scn, scn->kdtree, 0)) { 5.126 + 5.127 + ldir.normalize(); 5.128 + 5.129 + Vector3 vdir = -raydir / RAY_MAG; 5.130 + Vector3 vref = reflect(vdir, norm); 5.131 + 5.132 + float ndotl = dot(ldir, norm); 5.133 + float diff = MAX(ndotl, 0.0f); 5.134 + 5.135 + dcol[0] += mat->kd[0] * diff; 5.136 + dcol[1] += mat->kd[1] * diff; 5.137 + dcol[2] += mat->kd[2] * diff; 5.138 + 5.139 + float ldotvr = dot(ldir, vref); 5.140 + float spec = pow(MAX(ldotvr, 0.0f), mat->spow); 5.141 + 5.142 + scol[0] += mat->ks[0] * spec; 5.143 + scol[1] += mat->ks[1] * spec; 5.144 + scol[2] += mat->ks[2] * spec; 5.145 + } 5.146 + } 5.147 + 5.148 + float refl_color[3]; 5.149 + refl_color[0] = mat->ks[0] * mat->kr; 5.150 + refl_color[1] = mat->ks[1] * mat->kr; 5.151 + refl_color[2] = mat->ks[2] * mat->kr; 5.152 + 5.153 + energy *= (refl_color[0] + refl_color[1] + refl_color[2]) / 3.0; 5.154 + if(iter >= 0 && energy > MIN_ENERGY) { 5.155 + Vector3 rdir = reflect(-raydir, norm); 5.156 + 5.157 + Ray refl; 5.158 + refl.origin[0] = sp.pos.x; 5.159 + refl.origin[1] = sp.pos.y; 5.160 + refl.origin[2] = sp.pos.z; 5.161 + refl.dir[0] = rdir.x; 5.162 + refl.dir[1] = rdir.y; 5.163 + refl.dir[2] = rdir.z; 5.164 + 5.165 + float rcol[3]; 5.166 + trace_ray(rcol, refl, iter - 1, energy); 5.167 + scol[0] += rcol[0] * mat->ks[0] * mat->kr; 5.168 + scol[1] += rcol[1] * mat->ks[1] * mat->kr; 5.169 + scol[2] += rcol[2] * mat->ks[2] * mat->kr; 5.170 + } 5.171 + 5.172 + pixel[0] = dcol[0] + scol[0]; 5.173 + pixel[1] = dcol[1] + scol[1]; 5.174 + pixel[2] = dcol[2] + scol[2]; 5.175 +} 5.176 + 5.177 +static bool find_intersection(const Ray &ray, const Scene *scn, const KDNode *kd, SurfPoint *spret) 5.178 +{ 5.179 + if(!ray_aabb_test(ray, kd->aabb)) { 5.180 + return false; 5.181 + } 5.182 + 5.183 + SurfPoint sp, sptmp; 5.184 + if(!spret) { 5.185 + spret = &sptmp; 5.186 + } 5.187 + 5.188 + spret->t = RAY_MAG; 5.189 + spret->face = 0; 5.190 + 5.191 + if(kd->left) { 5.192 + assert(kd->right); 5.193 + 5.194 + bool found = find_intersection(ray, scn, kd->left, spret); 5.195 + if(find_intersection(ray, scn, kd->right, &sp)) { 5.196 + if(!found || sp.t < spret->t) { 5.197 + *spret = sp; 5.198 + } 5.199 + found = true; 5.200 + } 5.201 + return found; 5.202 + } 5.203 + 5.204 + const Face *faces = scn->get_face_buffer(); 5.205 + 5.206 + for(size_t i=0; i<kd->face_idx.size(); i++) { 5.207 + if(ray_triangle_test(ray, faces + kd->face_idx[i], &sp) && sp.t < spret->t) { 5.208 + *spret = sp; 5.209 + } 5.210 + } 5.211 + return spret->face != 0; 5.212 +} 5.213 + 5.214 +static bool ray_aabb_test(const Ray &ray, const AABBox &aabb) 5.215 +{ 5.216 + if(ray.origin[0] >= aabb.min[0] && ray.origin[1] >= aabb.min[1] && ray.origin[2] >= aabb.min[2] && 5.217 + ray.origin[0] < aabb.max[0] && ray.origin[1] < aabb.max[1] && ray.origin[2] < aabb.max[2]) { 5.218 + return true; 5.219 + } 5.220 + 5.221 + float bbox[][3] = { 5.222 + {aabb.min[0], aabb.min[1], aabb.min[2]}, 5.223 + {aabb.max[0], aabb.max[1], aabb.max[2]} 5.224 + }; 5.225 + 5.226 + int xsign = (int)(ray.dir[0] < 0.0); 5.227 + float invdirx = 1.0 / ray.dir[0]; 5.228 + float tmin = (bbox[xsign][0] - ray.origin[0]) * invdirx; 5.229 + float tmax = (bbox[1 - xsign][0] - ray.origin[0]) * invdirx; 5.230 + 5.231 + int ysign = (int)(ray.dir[1] < 0.0); 5.232 + float invdiry = 1.0 / ray.dir[1]; 5.233 + float tymin = (bbox[ysign][1] - ray.origin[1]) * invdiry; 5.234 + float tymax = (bbox[1 - ysign][1] - ray.origin[1]) * invdiry; 5.235 + 5.236 + if(tmin > tymax || tymin > tmax) { 5.237 + return false; 5.238 + } 5.239 + 5.240 + if(tymin > tmin) tmin = tymin; 5.241 + if(tymax < tmax) tmax = tymax; 5.242 + 5.243 + int zsign = (int)(ray.dir[2] < 0.0); 5.244 + float invdirz = 1.0 / ray.dir[2]; 5.245 + float tzmin = (bbox[zsign][2] - ray.origin[2]) * invdirz; 5.246 + float tzmax = (bbox[1 - zsign][2] - ray.origin[2]) * invdirz; 5.247 + 5.248 + if(tmin > tzmax || tzmin > tmax) { 5.249 + return false; 5.250 + } 5.251 + 5.252 + return tmin < 1.0 && tmax > 0.0; 5.253 + 5.254 +} 5.255 + 5.256 +static bool ray_triangle_test(const Ray &ray, const Face *face, SurfPoint *sp) 5.257 +{ 5.258 + Vector3 origin = ray.origin; 5.259 + Vector3 dir = ray.dir; 5.260 + Vector3 norm = face->normal; 5.261 + 5.262 + float ndotdir = dot(dir, norm); 5.263 + 5.264 + if(fabs(ndotdir) <= EPSILON) { 5.265 + return false; 5.266 + } 5.267 + 5.268 + Vector3 pt = face->v[0].pos; 5.269 + Vector3 vec = pt - origin; 5.270 + 5.271 + float ndotvec = dot(norm, vec); 5.272 + float t = ndotvec / ndotdir; 5.273 + 5.274 + if(t < EPSILON || t > 1.0) { 5.275 + return false; 5.276 + } 5.277 + pt = origin + dir * t; 5.278 + 5.279 + 5.280 + Vector3 bc = calc_bary(pt, face, norm); 5.281 + float bc_sum = bc.x + bc.y + bc.z; 5.282 + 5.283 + if(bc_sum < 1.0 - EPSILON || bc_sum > 1.0 + EPSILON) { 5.284 + return false; 5.285 + } 5.286 + 5.287 + Vector3 n0(face->v[0].normal); 5.288 + Vector3 n1(face->v[1].normal); 5.289 + Vector3 n2(face->v[2].normal); 5.290 + 5.291 + sp->t = t; 5.292 + sp->pos = pt; 5.293 + sp->norm = n0 * bc.x + n1 * bc.y + n2 * bc.z; 5.294 + sp->norm.normalize(); 5.295 + sp->face = face; 5.296 + return true; 5.297 +} 5.298 + 5.299 +static Vector3 calc_bary(const Vector3 &pt, const Face *face, const Vector3 &norm) 5.300 +{ 5.301 + Vector3 bc(0.0f, 0.0f, 0.0f); 5.302 + 5.303 + Vector3 v1 = Vector3(face->v[1].pos) - Vector3(face->v[0].pos); 5.304 + Vector3 v2 = Vector3(face->v[2].pos) - Vector3(face->v[0].pos); 5.305 + Vector3 xv1v2 = cross(v1, v2); 5.306 + 5.307 + float area = fabs(dot(xv1v2, norm)) * 0.5; 5.308 + if(area < EPSILON) { 5.309 + return bc; 5.310 + } 5.311 + 5.312 + Vector3 pv0 = face->v[0].pos - pt; 5.313 + Vector3 pv1 = face->v[1].pos - pt; 5.314 + Vector3 pv2 = face->v[2].pos - pt; 5.315 + 5.316 + // calculate the area of each sub-triangle 5.317 + Vector3 x12 = cross(pv1, pv2); 5.318 + Vector3 x20 = cross(pv2, pv0); 5.319 + Vector3 x01 = cross(pv0, pv1); 5.320 + 5.321 + float a0 = fabs(dot(x12, norm)) * 0.5; 5.322 + float a1 = fabs(dot(x20, norm)) * 0.5; 5.323 + float a2 = fabs(dot(x01, norm)) * 0.5; 5.324 + 5.325 + bc.x = a0 / area; 5.326 + bc.y = a1 / area; 5.327 + bc.z = a2 / area; 5.328 + return bc; 5.329 + 5.330 +} 5.331 + 5.332 +static void transform(float *res, const float *v, const float *xform) 5.333 +{ 5.334 + float tmp[3]; 5.335 + tmp[0] = v[0] * xform[0] + v[1] * xform[4] + v[2] * xform[8] + xform[12]; 5.336 + tmp[1] = v[0] * xform[1] + v[1] * xform[5] + v[2] * xform[9] + xform[13]; 5.337 + tmp[2] = v[0] * xform[2] + v[1] * xform[6] + v[2] * xform[10] + xform[14]; 5.338 + memcpy(res, tmp, sizeof tmp); 5.339 +} 5.340 + 5.341 +static void transform_ray(Ray *ray, const float *xform, const float *invtrans_xform) 5.342 +{ 5.343 + transform(ray->origin, ray->origin, xform); 5.344 + transform(ray->dir, ray->dir, invtrans_xform); 5.345 +}
6.1 --- a/src/rt.cc Sun Sep 05 16:43:55 2010 +0100 6.2 +++ b/src/rt.cc Fri Sep 10 16:47:00 2010 +0100 6.3 @@ -24,22 +24,6 @@ 6.4 NUM_KERNEL_ARGS 6.5 }; 6.6 6.7 -struct RendInfo { 6.8 - float ambient[4]; 6.9 - int xsz, ysz; 6.10 - int num_faces, num_lights; 6.11 - int max_iter; 6.12 - int cast_shadows; 6.13 -}; 6.14 - 6.15 -struct Ray { 6.16 - float origin[4], dir[4]; 6.17 -}; 6.18 - 6.19 -struct Light { 6.20 - float pos[4], color[4]; 6.21 -}; 6.22 - 6.23 static void update_render_info(); 6.24 static Ray get_primary_ray(int x, int y, int w, int h, float vfov_deg); 6.25 static float *create_kdimage(const KDNodeGPU *kdtree, int num_nodes, int *xsz_ret, int *ysz_ret); 6.26 @@ -49,10 +33,6 @@ 6.27 static CLProgram *prog; 6.28 static int global_size; 6.29 6.30 -static Light lightlist[] = { 6.31 - {{-8, 15, 18, 0}, {1, 1, 1, 1}} 6.32 -}; 6.33 - 6.34 6.35 static RendInfo rinf; 6.36 static int saved_iter_val; 6.37 @@ -72,7 +52,7 @@ 6.38 rinf.xsz = xsz; 6.39 rinf.ysz = ysz; 6.40 rinf.num_faces = scn->get_num_faces(); 6.41 - rinf.num_lights = sizeof lightlist / sizeof *lightlist; 6.42 + rinf.num_lights = scn->get_num_lights(); 6.43 rinf.max_iter = saved_iter_val = 6; 6.44 rinf.cast_shadows = true; 6.45 6.46 @@ -84,10 +64,11 @@ 6.47 prim_rays[i * xsz + j] = get_primary_ray(j, i, xsz, ysz, 45.0); 6.48 } 6.49 } 6.50 + dbg_set_primary_rays(prim_rays); // give them to the debug renderer 6.51 6.52 /* setup opencl */ 6.53 prog = new CLProgram("render"); 6.54 - if(!prog->load("rt.cl")) { 6.55 + if(!prog->load("src/rt.cl")) { 6.56 return false; 6.57 } 6.58 6.59 @@ -114,7 +95,7 @@ 6.60 prog->set_arg_buffer(KARG_RENDER_INFO, ARG_RD, sizeof rinf, &rinf); 6.61 prog->set_arg_buffer(KARG_FACES, ARG_RD, rinf.num_faces * sizeof(Face), faces); 6.62 prog->set_arg_buffer(KARG_MATLIB, ARG_RD, scn->get_num_materials() * sizeof(Material), scn->get_materials()); 6.63 - prog->set_arg_buffer(KARG_LIGHTS, ARG_RD, sizeof lightlist, lightlist); 6.64 + prog->set_arg_buffer(KARG_LIGHTS, ARG_RD, scn->get_num_lights() * sizeof(Light), scn->get_lights()); 6.65 prog->set_arg_buffer(KARG_PRIM_RAYS, ARG_RD, xsz * ysz * sizeof *prim_rays, prim_rays); 6.66 prog->set_arg_buffer(KARG_XFORM, ARG_RD, 16 * sizeof(float)); 6.67 prog->set_arg_buffer(KARG_INVTRANS_XFORM, ARG_RD, 16 * sizeof(float)); 6.68 @@ -133,9 +114,12 @@ 6.69 return false; 6.70 } 6.71 6.72 - delete [] prim_rays; 6.73 + //delete [] prim_rays; now dbg_renderer handles them 6.74 6.75 global_size = xsz * ysz; 6.76 + 6.77 + 6.78 + init_dbg_renderer(xsz, ysz, scn, tex); 6.79 return true; 6.80 } 6.81 6.82 @@ -143,7 +127,11 @@ 6.83 { 6.84 delete prog; 6.85 6.86 - printf("rendertime mean: %ld msec\n", timing_sample_sum / num_timing_samples); 6.87 + destroy_dbg_renderer(); 6.88 + 6.89 + if(num_timing_samples) { 6.90 + printf("rendertime mean: %ld msec\n", timing_sample_sum / num_timing_samples); 6.91 + } 6.92 } 6.93 6.94 bool render() 6.95 @@ -201,75 +189,6 @@ 6.96 return true; 6.97 } 6.98 6.99 -#define MIN(a, b) ((a) < (b) ? (a) : (b)) 6.100 -static void dbg_set_gl_material(Material *mat) 6.101 -{ 6.102 - static Material def_mat = {{0.7, 0.7, 0.7, 1}, {0, 0, 0, 0}, 0, 0, 0}; 6.103 - 6.104 - if(!mat) mat = &def_mat; 6.105 - 6.106 - glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat->kd); 6.107 - glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat->ks); 6.108 - glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, MIN(mat->spow, 128.0f)); 6.109 -} 6.110 - 6.111 -void dbg_render_gl(Scene *scn, bool show_tree, bool show_obj) 6.112 -{ 6.113 - glPushAttrib(GL_ENABLE_BIT | GL_TRANSFORM_BIT | GL_LIGHTING_BIT); 6.114 - 6.115 - for(int i=0; i<rinf.num_lights; i++) { 6.116 - float lpos[4]; 6.117 - 6.118 - memcpy(lpos, lightlist[i].pos, sizeof lpos); 6.119 - lpos[3] = 1.0; 6.120 - 6.121 - glLightfv(GL_LIGHT0 + i, GL_POSITION, lpos); 6.122 - glLightfv(GL_LIGHT0 + i, GL_DIFFUSE, lightlist[i].color); 6.123 - glEnable(GL_LIGHT0 + i); 6.124 - } 6.125 - 6.126 - glDisable(GL_TEXTURE_2D); 6.127 - glEnable(GL_DEPTH_TEST); 6.128 - glEnable(GL_LIGHTING); 6.129 - 6.130 - glMatrixMode(GL_PROJECTION); 6.131 - glPushMatrix(); 6.132 - glLoadIdentity(); 6.133 - gluPerspective(45.0, (float)rinf.xsz / (float)rinf.ysz, 0.5, 1000.0); 6.134 - 6.135 - if(show_obj) { 6.136 - Material *materials = scn->get_materials(); 6.137 - 6.138 - int num_faces = scn->get_num_faces(); 6.139 - int cur_mat = -1; 6.140 - 6.141 - for(int i=0; i<num_faces; i++) { 6.142 - if(faces[i].matid != cur_mat) { 6.143 - if(cur_mat != -1) { 6.144 - glEnd(); 6.145 - } 6.146 - dbg_set_gl_material(materials ? materials + faces[i].matid : 0); 6.147 - cur_mat = faces[i].matid; 6.148 - glBegin(GL_TRIANGLES); 6.149 - } 6.150 - 6.151 - for(int j=0; j<3; j++) { 6.152 - glNormal3fv(faces[i].v[j].normal); 6.153 - glVertex3fv(faces[i].v[j].pos); 6.154 - } 6.155 - } 6.156 - glEnd(); 6.157 - } 6.158 - 6.159 - if(show_tree) { 6.160 - scn->draw_kdtree(); 6.161 - } 6.162 - 6.163 - glPopMatrix(); 6.164 - glPopAttrib(); 6.165 - 6.166 - assert(glGetError() == GL_NO_ERROR); 6.167 -} 6.168 6.169 void set_xform(float *matrix, float *invtrans) 6.170 { 6.171 @@ -286,6 +205,12 @@ 6.172 unmap_mem_buffer(mbuf_invtrans); 6.173 } 6.174 6.175 + 6.176 +const RendInfo *get_render_info() 6.177 +{ 6.178 + return &rinf; 6.179 +} 6.180 + 6.181 void set_render_option(int opt, bool val) 6.182 { 6.183 switch(opt) { 6.184 @@ -403,6 +328,8 @@ 6.185 return ray; 6.186 } 6.187 6.188 +#define MIN(a, b) ((a) < (b) ? (a) : (b)) 6.189 + 6.190 static float *create_kdimage(const KDNodeGPU *kdtree, int num_nodes, int *xsz_ret, int *ysz_ret) 6.191 { 6.192 int ysz = MIN(num_nodes, KDIMG_MAX_HEIGHT);
7.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 7.2 +++ b/src/rt.cl Fri Sep 10 16:47:00 2010 +0100 7.3 @@ -0,0 +1,393 @@ 7.4 +/* vim: set ft=opencl:ts=4:sw=4 */ 7.5 +#include "common.h" 7.6 + 7.7 +struct RendInfo { 7.8 + float4 ambient; 7.9 + int xsz, ysz; 7.10 + int num_faces, num_lights; 7.11 + int max_iter; 7.12 + int cast_shadows; 7.13 +}; 7.14 + 7.15 +struct Vertex { 7.16 + float4 pos; 7.17 + float4 normal; 7.18 + float4 tex; 7.19 + float4 padding; 7.20 +}; 7.21 + 7.22 +struct Face { 7.23 + struct Vertex v[3]; 7.24 + float4 normal; 7.25 + int matid; 7.26 + int padding[3]; 7.27 +}; 7.28 + 7.29 +struct Material { 7.30 + float4 kd, ks; 7.31 + float kr, kt; 7.32 + float spow; 7.33 + float padding; 7.34 +}; 7.35 + 7.36 +struct Light { 7.37 + float4 pos, color; 7.38 +}; 7.39 + 7.40 +struct Ray { 7.41 + float4 origin, dir; 7.42 +}; 7.43 + 7.44 +struct SurfPoint { 7.45 + float t; 7.46 + float4 pos, norm, dbg; 7.47 + global const struct Face *obj; 7.48 + struct Material mat; 7.49 +}; 7.50 + 7.51 +struct Scene { 7.52 + float4 ambient; 7.53 + global const struct Face *faces; 7.54 + int num_faces; 7.55 + global const struct Light *lights; 7.56 + int num_lights; 7.57 + global const struct Material *matlib; 7.58 + //global const struct KDNode *kdtree; 7.59 + bool cast_shadows; 7.60 +}; 7.61 + 7.62 +struct AABBox { 7.63 + float4 min, max; 7.64 +}; 7.65 + 7.66 +struct KDNode { 7.67 + struct AABBox aabb; 7.68 + int face_idx[MAX_NODE_FACES]; 7.69 + int num_faces; 7.70 + int left, right; 7.71 + int padding; 7.72 +}; 7.73 + 7.74 +#define MIN_ENERGY 0.001 7.75 + 7.76 +float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg); 7.77 +bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *sp, read_only image2d_t kdimg); 7.78 +bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp); 7.79 +bool intersect_aabb(struct Ray ray, struct AABBox aabb); 7.80 + 7.81 +float4 reflect(float4 v, float4 n); 7.82 +float4 transform(float4 v, global const float *xform); 7.83 +void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans); 7.84 +float4 calc_bary(float4 pt, global const struct Face *face, float4 norm); 7.85 +float mean(float4 v); 7.86 + 7.87 +void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg); 7.88 + 7.89 + 7.90 +kernel void render(write_only image2d_t fb, 7.91 + global const struct RendInfo *rinf, 7.92 + global const struct Face *faces, 7.93 + global const struct Material *matlib, 7.94 + global const struct Light *lights, 7.95 + global const struct Ray *primrays, 7.96 + global const float *xform, 7.97 + global const float *invtrans, 7.98 + //global const struct KDNode *kdtree 7.99 + read_only image2d_t kdtree_img) 7.100 +{ 7.101 + int idx = get_global_id(0); 7.102 + 7.103 + struct Scene scn; 7.104 + scn.ambient = rinf->ambient; 7.105 + scn.faces = faces; 7.106 + scn.num_faces = rinf->num_faces; 7.107 + scn.lights = lights; 7.108 + scn.num_lights = rinf->num_lights; 7.109 + scn.matlib = matlib; 7.110 + scn.cast_shadows = rinf->cast_shadows; 7.111 + 7.112 + struct Ray ray = primrays[idx]; 7.113 + transform_ray(&ray, xform, invtrans); 7.114 + 7.115 + float4 pixel = (float4)(0, 0, 0, 0); 7.116 + float4 energy = (float4)(1.0, 1.0, 1.0, 0.0); 7.117 + int iter = 0; 7.118 + 7.119 + while(iter++ <= rinf->max_iter && mean(energy) > MIN_ENERGY) { 7.120 + struct SurfPoint sp; 7.121 + if(find_intersection(ray, &scn, &sp, kdtree_img)) { 7.122 + pixel += shade(ray, &scn, &sp, kdtree_img) * energy; 7.123 + 7.124 + float4 refl_col = sp.mat.ks * sp.mat.kr; 7.125 + 7.126 + ray.origin = sp.pos; 7.127 + ray.dir = reflect(-ray.dir, sp.norm); 7.128 + 7.129 + energy *= refl_col; 7.130 + } else { 7.131 + energy = (float4)(0.0, 0.0, 0.0, 0.0); 7.132 + } 7.133 + } 7.134 + 7.135 + int2 coord; 7.136 + coord.x = idx % rinf->xsz; 7.137 + coord.y = idx / rinf->xsz; 7.138 + 7.139 + write_imagef(fb, coord, pixel); 7.140 +} 7.141 + 7.142 +float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg) 7.143 +{ 7.144 + float4 norm = sp->norm; 7.145 + //bool entering = true; 7.146 + 7.147 + if(dot(ray.dir, norm) >= 0.0) { 7.148 + norm = -norm; 7.149 + //entering = false; 7.150 + } 7.151 + 7.152 + float4 dcol = scn->ambient * sp->mat.kd; 7.153 + float4 scol = (float4)(0, 0, 0, 0); 7.154 + 7.155 + for(int i=0; i<scn->num_lights; i++) { 7.156 + float4 ldir = scn->lights[i].pos - sp->pos; 7.157 + 7.158 + struct Ray shadowray; 7.159 + shadowray.origin = sp->pos; 7.160 + shadowray.dir = ldir; 7.161 + 7.162 + if(!scn->cast_shadows || !find_intersection(shadowray, scn, 0, kdimg)) { 7.163 + ldir = normalize(ldir); 7.164 + float4 vdir = -ray.dir; 7.165 + vdir.x = native_divide(vdir.x, RAY_MAG); 7.166 + vdir.y = native_divide(vdir.y, RAY_MAG); 7.167 + vdir.z = native_divide(vdir.z, RAY_MAG); 7.168 + float4 vref = reflect(vdir, norm); 7.169 + 7.170 + float diff = fmax(dot(ldir, norm), 0.0f); 7.171 + dcol += sp->mat.kd /* scn->lights[i].color*/ * diff; 7.172 + 7.173 + float spec = native_powr(fmax(dot(ldir, vref), 0.0f), sp->mat.spow); 7.174 + scol += sp->mat.ks /* scn->lights[i].color*/ * spec; 7.175 + } 7.176 + } 7.177 + 7.178 + return dcol + scol; 7.179 +} 7.180 + 7.181 +#define STACK_SIZE MAX_TREE_DEPTH 7.182 +bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *spres, read_only image2d_t kdimg) 7.183 +{ 7.184 + struct SurfPoint sp0; 7.185 + sp0.t = 1.0; 7.186 + sp0.obj = 0; 7.187 + 7.188 + int idxstack[STACK_SIZE]; 7.189 + int top = 0; // points after the topmost element of the stack 7.190 + idxstack[top++] = 0; // root at tree[0] 7.191 + 7.192 + while(top > 0) { 7.193 + int idx = idxstack[--top]; // remove this index from the stack and process it 7.194 + 7.195 + struct KDNode node; 7.196 + read_kdnode(idx, &node, kdimg); 7.197 + 7.198 + if(intersect_aabb(ray, node.aabb)) { 7.199 + if(node.left == -1) { 7.200 + // leaf node... check each face in turn and update the nearest intersection as needed 7.201 + for(int i=0; i<node.num_faces; i++) { 7.202 + struct SurfPoint spt; 7.203 + int fidx = node.face_idx[i]; 7.204 + 7.205 + if(intersect(ray, scn->faces + fidx, &spt) && spt.t < sp0.t) { 7.206 + sp0 = spt; 7.207 + } 7.208 + } 7.209 + } else { 7.210 + // internal node... recurse to the children 7.211 + idxstack[top++] = node.left; 7.212 + idxstack[top++] = node.right; 7.213 + } 7.214 + } 7.215 + } 7.216 + 7.217 + if(!sp0.obj) { 7.218 + return false; 7.219 + } 7.220 + 7.221 + if(spres) { 7.222 + *spres = sp0; 7.223 + spres->mat = scn->matlib[sp0.obj->matid]; 7.224 + } 7.225 + return true; 7.226 +} 7.227 + 7.228 +bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp) 7.229 +{ 7.230 + float4 origin = ray.origin; 7.231 + float4 dir = ray.dir; 7.232 + float4 norm = face->normal; 7.233 + 7.234 + float ndotdir = dot(dir, norm); 7.235 + 7.236 + if(fabs(ndotdir) <= EPSILON) { 7.237 + return false; 7.238 + } 7.239 + 7.240 + float4 pt = face->v[0].pos; 7.241 + float4 vec = pt - origin; 7.242 + 7.243 + float ndotvec = dot(norm, vec); 7.244 + float t = native_divide(ndotvec, ndotdir); 7.245 + 7.246 + if(t < EPSILON || t > 1.0) { 7.247 + return false; 7.248 + } 7.249 + pt = origin + dir * t; 7.250 + 7.251 + 7.252 + float4 bc = calc_bary(pt, face, norm); 7.253 + float bc_sum = bc.x + bc.y + bc.z; 7.254 + 7.255 + if(bc_sum < 1.0 - EPSILON || bc_sum > 1.0 + EPSILON) { 7.256 + return false; 7.257 + } 7.258 + 7.259 + sp->t = t; 7.260 + sp->pos = pt; 7.261 + sp->norm = normalize(face->v[0].normal * bc.x + face->v[1].normal * bc.y + face->v[2].normal * bc.z); 7.262 + sp->obj = face; 7.263 + sp->dbg = bc; 7.264 + return true; 7.265 +} 7.266 + 7.267 +bool intersect_aabb(struct Ray ray, struct AABBox aabb) 7.268 +{ 7.269 + if(ray.origin.x >= aabb.min.x && ray.origin.y >= aabb.min.y && ray.origin.z >= aabb.min.z && 7.270 + ray.origin.x < aabb.max.x && ray.origin.y < aabb.max.y && ray.origin.z < aabb.max.z) { 7.271 + return true; 7.272 + } 7.273 + 7.274 + float4 bbox[2] = { 7.275 + aabb.min.x, aabb.min.y, aabb.min.z, 0, 7.276 + aabb.max.x, aabb.max.y, aabb.max.z, 0 7.277 + }; 7.278 + 7.279 + int xsign = (int)(ray.dir.x < 0.0); 7.280 + float invdirx = native_recip(ray.dir.x); 7.281 + float tmin = (bbox[xsign].x - ray.origin.x) * invdirx; 7.282 + float tmax = (bbox[1 - xsign].x - ray.origin.x) * invdirx; 7.283 + 7.284 + int ysign = (int)(ray.dir.y < 0.0); 7.285 + float invdiry = native_recip(ray.dir.y); 7.286 + float tymin = (bbox[ysign].y - ray.origin.y) * invdiry; 7.287 + float tymax = (bbox[1 - ysign].y - ray.origin.y) * invdiry; 7.288 + 7.289 + if(tmin > tymax || tymin > tmax) { 7.290 + return false; 7.291 + } 7.292 + 7.293 + if(tymin > tmin) tmin = tymin; 7.294 + if(tymax < tmax) tmax = tymax; 7.295 + 7.296 + int zsign = (int)(ray.dir.z < 0.0); 7.297 + float invdirz = native_recip(ray.dir.z); 7.298 + float tzmin = (bbox[zsign].z - ray.origin.z) * invdirz; 7.299 + float tzmax = (bbox[1 - zsign].z - ray.origin.z) * invdirz; 7.300 + 7.301 + if(tmin > tzmax || tzmin > tmax) { 7.302 + return false; 7.303 + } 7.304 + 7.305 + return tmin < 1.0 && tmax > 0.0; 7.306 +} 7.307 + 7.308 +float4 reflect(float4 v, float4 n) 7.309 +{ 7.310 + return 2.0f * dot(v, n) * n - v; 7.311 +} 7.312 + 7.313 +float4 transform(float4 v, global const float *xform) 7.314 +{ 7.315 + float4 res; 7.316 + res.x = v.x * xform[0] + v.y * xform[4] + v.z * xform[8] + xform[12]; 7.317 + res.y = v.x * xform[1] + v.y * xform[5] + v.z * xform[9] + xform[13]; 7.318 + res.z = v.x * xform[2] + v.y * xform[6] + v.z * xform[10] + xform[14]; 7.319 + res.w = 0.0; 7.320 + return res; 7.321 +} 7.322 + 7.323 +void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans) 7.324 +{ 7.325 + ray->origin = transform(ray->origin, xform); 7.326 + ray->dir = transform(ray->dir, invtrans); 7.327 +} 7.328 + 7.329 +float4 calc_bary(float4 pt, global const struct Face *face, float4 norm) 7.330 +{ 7.331 + float4 bc = (float4)(0, 0, 0, 0); 7.332 + 7.333 + // calculate area of the whole triangle 7.334 + float4 v1 = face->v[1].pos - face->v[0].pos; 7.335 + float4 v2 = face->v[2].pos - face->v[0].pos; 7.336 + float4 xv1v2 = cross(v1, v2); 7.337 + 7.338 + float area = fabs(dot(xv1v2, norm)) * 0.5; 7.339 + if(area < EPSILON) { 7.340 + return bc; 7.341 + } 7.342 + 7.343 + float4 pv0 = face->v[0].pos - pt; 7.344 + float4 pv1 = face->v[1].pos - pt; 7.345 + float4 pv2 = face->v[2].pos - pt; 7.346 + 7.347 + // calculate the area of each sub-triangle 7.348 + float4 x12 = cross(pv1, pv2); 7.349 + float4 x20 = cross(pv2, pv0); 7.350 + float4 x01 = cross(pv0, pv1); 7.351 + 7.352 + float a0 = fabs(dot(x12, norm)) * 0.5; 7.353 + float a1 = fabs(dot(x20, norm)) * 0.5; 7.354 + float a2 = fabs(dot(x01, norm)) * 0.5; 7.355 + 7.356 + bc.x = native_divide(a0, area); 7.357 + bc.y = native_divide(a1, area); 7.358 + bc.z = native_divide(a2, area); 7.359 + return bc; 7.360 +} 7.361 + 7.362 +float mean(float4 v) 7.363 +{ 7.364 + return native_divide(v.x + v.y + v.z, 3.0); 7.365 +} 7.366 + 7.367 + 7.368 +const sampler_t kdsampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_NONE | CLK_FILTER_NEAREST; 7.369 + 7.370 +// read a KD-tree node from a texture scanline 7.371 +void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg) 7.372 +{ 7.373 + int startx = KDIMG_NODE_WIDTH * (idx / KDIMG_MAX_HEIGHT); 7.374 + 7.375 + int2 tc; 7.376 + tc.x = startx; 7.377 + tc.y = idx % KDIMG_MAX_HEIGHT; 7.378 + 7.379 + node->aabb.min = read_imagef(kdimg, kdsampler, tc); tc.x++; 7.380 + node->aabb.max = read_imagef(kdimg, kdsampler, tc); 7.381 + 7.382 + tc.x = startx + 2 + MAX_NODE_FACES / 4; 7.383 + float4 pix = read_imagef(kdimg, kdsampler, tc); 7.384 + node->num_faces = (int)pix.x; 7.385 + node->left = (int)pix.y; 7.386 + node->right = (int)pix.z; 7.387 + 7.388 + tc.x = startx + 2; 7.389 + for(int i=0; i<node->num_faces; i+=4) { 7.390 + float4 pix = read_imagef(kdimg, kdsampler, tc); tc.x++; 7.391 + node->face_idx[i] = (int)pix.x; 7.392 + node->face_idx[i + 1] = (int)pix.y; 7.393 + node->face_idx[i + 2] = (int)pix.z; 7.394 + node->face_idx[i + 3] = (int)pix.w; 7.395 + } 7.396 +}
8.1 --- a/src/rt.h Sun Sep 05 16:43:55 2010 +0100 8.2 +++ b/src/rt.h Fri Sep 10 16:47:00 2010 +0100 8.3 @@ -11,11 +11,26 @@ 8.4 NUM_RENDER_OPTIONS 8.5 }; 8.6 8.7 +struct RendInfo { 8.8 + float ambient[4]; 8.9 + int xsz, ysz; 8.10 + int num_faces, num_lights; 8.11 + int max_iter; 8.12 + int cast_shadows; 8.13 +}; 8.14 + 8.15 +struct Ray { 8.16 + float origin[4], dir[4]; 8.17 +}; 8.18 + 8.19 + 8.20 bool init_renderer(int xsz, int ysz, Scene *scn, unsigned int tex); 8.21 void destroy_renderer(); 8.22 bool render(); 8.23 void set_xform(float *matrix, float *invtrans); 8.24 8.25 +const RendInfo *get_render_info(); 8.26 + 8.27 void set_render_option(int opt, bool val); 8.28 void set_render_option(int opt, int val); 8.29 void set_render_option(int opt, float val); 8.30 @@ -24,6 +39,13 @@ 8.31 int get_render_option_int(int opt); 8.32 float get_render_option_float(int opt); 8.33 8.34 +// raytrace in the CPU 8.35 +bool init_dbg_renderer(int xsz, int ysz, Scene *scn, unsigned int texid); 8.36 +void destroy_dbg_renderer(); 8.37 +void dbg_set_primary_rays(const Ray *rays); 8.38 +void dbg_render(const float *xform, const float *invtrans_xform, int num_threads = -1); 8.39 + 8.40 +// visualize the scene using OpenGL 8.41 void dbg_render_gl(Scene *scn, bool show_tree = false, bool show_obj = true); 8.42 8.43 #endif /* RT_H_ */
9.1 --- a/src/scene.cc Sun Sep 05 16:43:55 2010 +0100 9.2 +++ b/src/scene.cc Fri Sep 10 16:47:00 2010 +0100 9.3 @@ -109,11 +109,27 @@ 9.4 return true; 9.5 } 9.6 9.7 +bool Scene::add_light(const Light <) 9.8 +{ 9.9 + try { 9.10 + lights.push_back(lt); 9.11 + } 9.12 + catch(...) { 9.13 + return false; 9.14 + } 9.15 + return true; 9.16 +} 9.17 + 9.18 int Scene::get_num_meshes() const 9.19 { 9.20 return (int)meshes.size(); 9.21 } 9.22 9.23 +int Scene::get_num_lights() const 9.24 +{ 9.25 + return (int)lights.size(); 9.26 +} 9.27 + 9.28 int Scene::get_num_faces() const 9.29 { 9.30 if(num_faces >= 0) { 9.31 @@ -137,6 +153,38 @@ 9.32 return kdtree_nodes(kdtree); 9.33 } 9.34 9.35 +Mesh **Scene::get_meshes() 9.36 +{ 9.37 + if(meshes.empty()) { 9.38 + return 0; 9.39 + } 9.40 + return &meshes[0]; 9.41 +} 9.42 + 9.43 +const Mesh * const *Scene::get_meshes() const 9.44 +{ 9.45 + if(meshes.empty()) { 9.46 + return 0; 9.47 + } 9.48 + return &meshes[0]; 9.49 +} 9.50 + 9.51 +Light *Scene::get_lights() 9.52 +{ 9.53 + if(lights.empty()) { 9.54 + return 0; 9.55 + } 9.56 + return &lights[0]; 9.57 +} 9.58 + 9.59 +const Light *Scene::get_lights() const 9.60 +{ 9.61 + if(lights.empty()) { 9.62 + return 0; 9.63 + } 9.64 + return &lights[0]; 9.65 +} 9.66 + 9.67 Material *Scene::get_materials() 9.68 { 9.69 if(matlib.empty()) {
10.1 --- a/src/scene.h Sun Sep 05 16:43:55 2010 +0100 10.2 +++ b/src/scene.h Fri Sep 10 16:47:00 2010 +0100 10.3 @@ -34,6 +34,10 @@ 10.4 int matid; 10.5 }; 10.6 10.7 +struct Light { 10.8 + float pos[4], color[4]; 10.9 +}; 10.10 + 10.11 class AABBox { 10.12 public: 10.13 float min[4], max[4]; 10.14 @@ -70,6 +74,7 @@ 10.15 10.16 public: 10.17 std::vector<Mesh*> meshes; 10.18 + std::vector<Light> lights; 10.19 std::vector<Material> matlib; 10.20 KDNode *kdtree; 10.21 10.22 @@ -77,11 +82,20 @@ 10.23 ~Scene(); 10.24 10.25 bool add_mesh(Mesh *m); 10.26 + bool add_light(const Light <); 10.27 + 10.28 int get_num_meshes() const; 10.29 + int get_num_lights() const; 10.30 int get_num_faces() const; 10.31 int get_num_materials() const; 10.32 int get_num_kdnodes() const; 10.33 10.34 + Mesh **get_meshes(); 10.35 + const Mesh * const *get_meshes() const; 10.36 + 10.37 + Light *get_lights(); 10.38 + const Light *get_lights() const; 10.39 + 10.40 Material *get_materials(); 10.41 const Material *get_materials() const; 10.42
11.1 --- a/src/vector.cc Sun Sep 05 16:43:55 2010 +0100 11.2 +++ b/src/vector.cc Fri Sep 10 16:47:00 2010 +0100 11.3 @@ -18,6 +18,13 @@ 11.4 this->z = z; 11.5 } 11.6 11.7 +Vector3::Vector3(const float *arr) 11.8 +{ 11.9 + x = arr[0]; 11.10 + y = arr[1]; 11.11 + z = arr[2]; 11.12 +} 11.13 + 11.14 void Vector3::normalize() 11.15 { 11.16 float len = sqrt(x * x + y * y + z * z);
12.1 --- a/src/vector.h Sun Sep 05 16:43:55 2010 +0100 12.2 +++ b/src/vector.h Fri Sep 10 16:47:00 2010 +0100 12.3 @@ -15,6 +15,7 @@ 12.4 12.5 Vector3(); 12.6 Vector3(float x, float y, float z); 12.7 + Vector3(const float *arr); 12.8 12.9 void normalize(); 12.10 inline float length(); 12.11 @@ -28,10 +29,13 @@ 12.12 12.13 inline Vector3 operator -(const Vector3 &vec); 12.14 inline Vector3 operator *(const Vector3 &vec, float s); 12.15 +inline Vector3 operator /(const Vector3 &vec, float s); 12.16 12.17 inline float dot(const Vector3 &a, const Vector3 &b); 12.18 inline Vector3 cross(const Vector3 &a, const Vector3 &b); 12.19 12.20 +inline Vector3 reflect(const Vector3 &v, const Vector3 &n); 12.21 + 12.22 #include "vector.inl" 12.23 12.24 #endif /* VECTOR_H_ */
13.1 --- a/src/vector.inl Sun Sep 05 16:43:55 2010 +0100 13.2 +++ b/src/vector.inl Fri Sep 10 16:47:00 2010 +0100 13.3 @@ -41,6 +41,10 @@ 13.4 return Vector3(vec.x * s, vec.y * s, vec.z * s); 13.5 } 13.6 13.7 +inline Vector3 operator /(const Vector3 &vec, float s) 13.8 +{ 13.9 + return Vector3(vec.x / s, vec.y / s, vec.z / s); 13.10 +} 13.11 13.12 inline float dot(const Vector3 &a, const Vector3 &b) 13.13 { 13.14 @@ -51,3 +55,8 @@ 13.15 { 13.16 return Vector3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); 13.17 } 13.18 + 13.19 +inline Vector3 reflect(const Vector3 &v, const Vector3 &n) 13.20 +{ 13.21 + return n * (2.0 * dot(v, n)) - v; 13.22 +}