cleaned Texture interface
new C++ demo: textures
slightly adjusted SAH for kd-tree
slightly optimized kd-tree building -- moved termination cond. so it's tested before recursion
minor sphere intersection optimization
/* * Pyrit Ray Tracer * file: raytracer.cc * * Radek Brich, 2006-2007 */#ifdef PTHREADS#include <pthread.h>#endif#include <stdio.h>#include <malloc.h>#include "raytracer.h"int pyrit_verbosity = 2;// Hammersley spherical point distribution// http://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoints.htmlVector3 Raytracer::SphereDistribute(int i, int n, Float extent, Vector3 &normal){ Float p, t, st, phi, phirad; int kk; t = 0; for (p=0.5, kk=i; kk; p*=0.5, kk>>=1) if (kk & 1) t += p; t = 1.0 + (t - 1.0)*extent; phi = (i + 0.5) / n; phirad = phi * 2.0 * M_PI; st = sqrt(1.0 - t*t); Float x, y, z, xx, yy, zz, q; x = st * cos(phirad); y = st * sin(phirad); z = t; // rotate against Y axis q = acos(normal.z); zz = z*cos(q) - x*sin(q); xx = z*sin(q) + x*cos(q); yy = y; // rotate against Z axis q = atan2f(normal.y, normal.x); x = xx*cos(q) - yy*sin(q); y = xx*sin(q) + yy*cos(q); z = zz; return Vector3(x, y, z);}// ---- tyto dve funkce budou v budouci verzi metody objektu PhongShader// calculate shader function// P is point of intersection, N normal in this pointColour PhongShader_ambient(Material &mat, Vector3 &P){ Colour col; if (mat.texture) col = mat.texture->evaluate(P); else col = mat.colour; // ambient return mat.ambient * col;}/* P is point of intersection, N normal in this point, R direction of reflected ray, V direction to the viewer*/Colour PhongShader_calculate(Material &mat, Vector3 &P, Vector3 &N, Vector3 &R, Vector3 &V, Light &light){ Colour I = Colour(); Vector3 L = light.pos - P; L.normalize(); Float L_dot_N = dot(L, N); Float R_dot_V = dot(R, V); Colour col; if (mat.texture) col = mat.texture->evaluate(P); else col = mat.colour; // diffuse I = mat.diffuse * col * light.colour * L_dot_N; // specular if (R_dot_V > 0) I += mat.specular * light.colour * powf(R_dot_V, mat.shininess); return I;}Colour Raytracer::raytrace(Ray &ray, int depth, Shape *origin_shape){ Float nearest_distance = Inf; Shape *nearest_shape = top->nearest_intersection(origin_shape, ray, nearest_distance); if (nearest_shape == NULL) { return bg_colour; } else { Colour col = Colour(); Vector3 P = ray.o + ray.dir * nearest_distance; // point of intersection Vector3 normal = nearest_shape->normal(P); bool from_inside = false; // make shapes double sided if (dot(normal, ray.dir) > 0.0) { normal = - normal; from_inside = true; } col = PhongShader_ambient(*nearest_shape->material, P); vector<Light*>::iterator light; for (light = lights.begin(); light != lights.end(); light++) { Vector3 jo, L = (*light)->pos - P; // direction vector to light L.normalize(); Float L_dot_N = dot(L, normal); if (L_dot_N > 0) { // test if this light is occluded (sharp shadows) if ((*light)->cast_shadows) { Ray shadow_ray = Ray(P, L); Float dist = FLT_MAX; if (top->nearest_intersection(nearest_shape, shadow_ray, dist)) continue; } // shading function Vector3 R = L - 2.0 * L_dot_N * normal; col += PhongShader_calculate(*nearest_shape->material, P, normal, R, ray.dir, **light); } } if (depth < max_depth) { Colour trans_col, refl_col; Float trans = nearest_shape->material->transmissivity; Float refl = nearest_shape->material->reflectivity; const Float cos_i = - dot(normal, ray.dir); // reflection if (refl > 0.01) { Vector3 newdir = ray.dir + 2.0 * cos_i * normal; Ray newray = Ray(P, newdir); refl_col = raytrace(newray, depth + 1, nearest_shape); } // refraction if (trans > 0.01) { Float n, n1, n2; if (from_inside) { n1 = nearest_shape->material->refract_index; n2 = 1.0; n = n1; } else { n1 = 1.0; n2 = nearest_shape->material->refract_index; n = 1.0 / n2; } const Float sin2_t = n*n * (1 - cos_i*cos_i); if (sin2_t >= 1.0) { // totally reflected refl += trans; trans = 0; } else { const Float cos_t = sqrtf(1 - sin2_t); const Float Rdiv = 1.0/(n1*cos_i + n2*cos_t); const Float Rper = (n1*cos_i - n2*cos_t)*Rdiv; const Float Rpar = (n2*cos_i - n1*cos_t)*Rdiv; const Float R = (Rper*Rper + Rpar*Rpar)/2; refl += R*trans; trans = (1-R)*trans; Vector3 newdir = n * ray.dir + (n*cos_i - cos_t) * normal; Ray newray = Ray(P, newdir); trans_col = raytrace(newray, depth + 1, nearest_shape); } } col = (1-refl-trans)*col + refl*refl_col + trans*trans_col; } // ambient occlusion if (!from_inside && ao_samples) { Float miss = 0; for (int i = 0; i < ao_samples; i++) { Vector3 dir = SphereDistribute(i, ao_samples, ao_angle, normal); Ray ao_ray = Ray(P, dir); Float dist = ao_distance; Shape *shape_in_way = top->nearest_intersection(nearest_shape, ao_ray, dist); if (shape_in_way == NULL) miss += 1.0; else miss += dist / ao_distance; } Float ao_intensity = miss / ao_samples; col = col * ao_intensity; } return col; }}static inline void samplepixel(Colour &c, Vector3 &dir, RenderrowData* d, const int &oversample){ if (oversample <= 0) { // no oversampling dir.normalize(); Ray ray(d->eye, dir); c = d->rt->raytrace(ray, 0, NULL); } else if (oversample <= 3) { // grid oversampling c = Colour(0,0,0); static const int gridsamples[] = {5,9,16}; static const Float osa5x[] = {0.0, -0.4, +0.4, +0.4, -0.4}; static const Float osa5y[] = {0.0, -0.4, -0.4, +0.4, +0.4}; static const Float osa9x[] = {-0.34, 0.00, +0.34, -0.34, 0.00, +0.34, -0.34, 0.00, +0.34}; static const Float osa9y[] = {-0.34, -0.34, -0.34, 0.00, 0.00, 0.00, +0.34, +0.34, +0.34}; static const Float osa16x[] = {-0.375, -0.125, +0.125, +0.375, -0.375, -0.125, +0.125, +0.375, -0.375, -0.125, +0.125, +0.375, -0.375, -0.125, +0.125, +0.375}; static const Float osa16y[] = {-0.375, -0.375, -0.375, -0.375, -0.125, -0.125, -0.125, -0.125, +0.125, +0.125, +0.125, +0.125, +0.375, +0.375, +0.375, +0.375}; static const Float *osaSx[] = {osa5x, osa9x, osa16x}; static const Float *osaSy[] = {osa5y, osa9y, osa16y}; const int samples = gridsamples[oversample-1]; const Float *osax = osaSx[oversample-1]; const Float *osay = osaSy[oversample-1]; for (int i = 0; i < samples; i++) { Vector3 tmpdir = dir + osax[i]*d->dx + osay[i]*d->dy; tmpdir.normalize(); Ray ray(d->eye, tmpdir); c += d->rt->raytrace(ray, 0, NULL); } c = c * (1.0/samples); } else { // stochastic oversampling // ...todo }}static void *renderrow(void *data){ RenderrowData *d = (RenderrowData*) data; const int subsample = d->rt->getSubsample(); const Float subsample2 = 1.0/(subsample*subsample); const int oversample = d->rt->getOversample(); const int ww = d->w*3; Vector3 dir = d->dfix; for (int x = 0; x < d->w; x += subsample) { // generate a ray from eye passing through this pixel if (subsample > 1) { Colour ic; // top-left dir.normalize(); Ray ray(d->eye, dir); Colour c1 = d->rt->raytrace(ray, 0, NULL); // top-right Vector3 tmpdir = dir + (subsample-1)*d->dx; tmpdir.normalize(); ray.dir = tmpdir; Colour c2 = d->rt->raytrace(ray, 0, NULL); // bottom right tmpdir += (subsample-1)*d->dy; tmpdir.normalize(); ray.dir = tmpdir; Colour c4 = d->rt->raytrace(ray, 0, NULL); // bottom left tmpdir = dir + (subsample-1)*d->dy; tmpdir.normalize(); ray.dir = tmpdir; Colour c3 = d->rt->raytrace(ray, 0, NULL); // are the colors similar? Float m = (c1-c2).mag2(); m = max(m, (c2-c3).mag2()); m = max(m, (c3-c4).mag2()); m = max(m, (c4-c1).mag2()); if (m < 0.001) { // interpolate Float *i = d->iter; for (int x = 0; x < subsample; x++) { for (int y = 0; y < subsample; y++) { ic = c1*(subsample-x)*(subsample-y)*subsample2 + c2*(x)*(subsample-y)*subsample2 + c3*(subsample-x)*(y)*subsample2 + c4*(x)*(y)*subsample2; *(i + ww*y) = ic.r; *(i + ww*y + 1) = ic.g; *(i + ww*y + 2) = ic.b; } i += 3; } d->iter = i; } else { // render all pixels Vector3 tmpdir = dir; if (oversample) { for (int x = 0; x < subsample; x++) { for (int y = 0; y < subsample; y++) { Vector3 tmp2dir = tmpdir + y*d->dy; samplepixel(ic, tmp2dir, d, oversample); *(d->iter + ww*y) = ic.r; *(d->iter + ww*y + 1) = ic.g; *(d->iter + ww*y + 2) = ic.b; } d->iter += 3; tmpdir += d->dx; } } else { /* this is so complex because it tries to reuse already computed corner pixels though, above code will also work for non-oversampling... */ // first column *(d->iter) = c1.r; *(d->iter + 1) = c1.g; *(d->iter + 2) = c1.b; for (int y = 1; y < subsample-1; y++) { Vector3 tmp2dir = tmpdir + y*d->dy; tmp2dir.normalize(); ray.dir = tmp2dir; ic = d->rt->raytrace(ray, 0, NULL); *(d->iter + ww*y) = ic.r; *(d->iter + ww*y + 1) = ic.g; *(d->iter + ww*y + 2) = ic.b; } *(d->iter + ww*(subsample-1)) = c3.r; *(d->iter + ww*(subsample-1) + 1) = c3.g; *(d->iter + ww*(subsample-1) + 2) = c3.b; d->iter += 3; tmpdir += d->dx; // middle for (int x = 1; x < subsample-1; x++) { for (int y = 0; y < subsample; y++) { Vector3 tmp2dir = tmpdir + y*d->dy; tmp2dir.normalize(); ray.dir = tmp2dir; ic = d->rt->raytrace(ray, 0, NULL); *(d->iter + ww*y) = ic.r; *(d->iter + ww*y + 1) = ic.g; *(d->iter + ww*y + 2) = ic.b; } d->iter += 3; tmpdir += d->dx; } // last column *(d->iter) = c2.r; *(d->iter + 1) = c2.g; *(d->iter + 2) = c2.b; for (int y = 1; y < subsample-1; y++) { Vector3 tmp2dir = tmpdir + y*d->dy; tmp2dir.normalize(); ray.dir = tmp2dir; ic = d->rt->raytrace(ray, 0, NULL); *(d->iter + ww*y) = ic.r; *(d->iter + ww*y + 1) = ic.g; *(d->iter + ww*y + 2) = ic.b; } *(d->iter + ww*(subsample-1)) = c4.r; *(d->iter + ww*(subsample-1) + 1) = c4.g; *(d->iter + ww*(subsample-1) + 2) = c4.b; d->iter += 3; } } } else // subsample <= 1 { Colour c; samplepixel(c, dir, d, oversample); // write color to buffer *d->iter++ = c.r; *d->iter++ = c.g; *d->iter++ = c.b; } dir += d->dx*subsample; }#ifdef PTHREADS pthread_exit((void *)d);#endif return (void *)d;}void Raytracer::render(int w, int h, Float *buffer){ if (!camera || !top || !buffer) return; RenderrowData *d; Float S = 0.5/w; Vector3 dfix = camera->u*(-w/2.0*S/camera->f) + camera->v*(h/2.0*S/camera->f) + camera->p; Vector3 dx = camera->u * (S/camera->f); Vector3 dy = camera->v * (-S/camera->f);#ifdef PTHREADS dbgmsg(1, "* pthreads enabled, using %d threads\n", num_threads); pthread_t threads[num_threads]; for (int t = 0; t < num_threads; t++) threads[t] = pthread_self(); int t = 0;#endif /* for each pixel... */ dbgmsg(1, "* raytracing...\n"); dbgmsg(2, "- row 0 ( 0%% done)"); for (int y = 0; y < h; y += subsample) { d = (RenderrowData*) malloc(sizeof(RenderrowData)); d->rt = this; d->w = w; d->eye = camera->eye; d->dfix = dfix; d->dx = dx; d->dy = dy; d->iter = buffer + y*3*w;#ifdef PTHREADS /* create new thread and increase 't' */ int rc = pthread_create(&threads[t++], NULL, renderrow, (void *)d); if (rc) { dbgmsg(0, "\nE pthread_create unsuccessful, return code was %d\n", rc); exit(1); } /* when 't' owerflows, reset it */ if (t >= num_threads) t = 0; /* wait for next thread in fifo queue, so the descriptor can be reused; this also limits number of running threads */ if (!pthread_equal(threads[t], pthread_self())) if (pthread_join(threads[t], (void**)&d) == 0) free(d);#else renderrow((void *)d); free(d);#endif dfix += dy*subsample; dbgmsg(2, "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%4d (%2d%% done)", y, y*100/(h-1)); } dbgmsg(2, "\n");#ifdef PTHREADS dbgmsg(2, "- waiting for threads to finish\n"); for (t = 0; t < num_threads; t++) if (pthread_join(threads[t], (void**)&d) == 0) free(d);#endif}void Raytracer::addlight(Light *light){ lights.push_back(light);}void Raytracer::ambientocclusion(int samples, Float distance, Float angle){ ao_samples = samples; ao_distance = distance; ao_angle = angle; if (ao_distance == 0) /* 0 ==> Inf */ ao_distance = FLT_MAX;}