more build script tuning
make all float constants single precision
solve many warnings from msvc and gcc with various -W... flags
add common.cc file for dbgmsg() function witch apparently cannot be inlined
fix python module building with msvc, add manifest file handling
remove forgotten RenderrowData class
add stanford models download script for windows (.bat)
/* * raytracer.cc: Raytracer class * * This file is part of Pyrit Ray Tracer. * * Copyright 2006, 2007, 2008 Radek Brich * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN * THE SOFTWARE. */#include <pthread.h>#include <stdio.h>#include <stdlib.h>#include <assert.h>#include "raytracer.h"// Hammersley spherical point distribution// http://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoints.htmlVector Raytracer::SphereDistribute(int i, int n, Float extent, const Vector &normal){ Float p, t, st, phi, phirad; int kk; t = 0; for (p=0.5f, kk=i; kk; p*=0.5f, kk>>=1) if (kk & 1) t += p; t = 1.0f + (t - 1.0f)*extent; phi = (i + 0.5f) / n; phirad = phi * 2.0f * PI; st = sqrt(1.0f - 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 Vector(x, y, z);}/* P is point of intersection, N normal in this point, R direction of reflected ray, V direction to the viewer*/Colour Raytracer::PhongShader(const Shape *shape, const Vector &P, const Vector &N, const Vector &V){ Colour col, acc; Material * const &mat = shape->material; if (mat->texture) col = mat->texture->evaluate(P); else col = mat->colour; // ambient acc = mat->ambient * col; vector<Light*>::iterator light; for (light = lights.begin(); light != lights.end(); light++) { const Vector L = normalize((*light)->pos - P); // direction vector to light const Float L_dot_N = dot(L, N); if (L_dot_N > 0) { // test if this light is occluded (sharp shadows) if ((*light)->cast_shadows) { const Ray shadow_ray = Ray(P, L); Float dist = FLT_MAX; if (top->nearest_intersection(shape, shadow_ray, dist)) continue; } const Vector R = L - 2.0f * L_dot_N * N; const Float R_dot_V = dot(R, V); // diffuse acc += mat->diffuse * col * (*light)->colour * L_dot_N; // specular if (R_dot_V > 0) acc += mat->specular * (*light)->colour * powf(R_dot_V, mat->shininess); } } return acc;}#ifndef NO_SIMDVectorPacket Raytracer::PhongShader_packet(const Shape* const* shapes, const VectorPacket &P, const VectorPacket &N, const VectorPacket &V){ VectorPacket acc, colour; union { mfloat4 ambient; float ambient_f[4]; }; union { mfloat4 diffuse; float diffuse_f[4]; }; union { mfloat4 specular; float specular_f[4]; }; union { mfloat4 shininess; float shininess_f[4]; }; for (int i = 0; i < 4; i++) if (shapes[i] == NULL) { ambient_f[i] = 0; diffuse_f[i] = 0; specular_f[i] = 0; shininess_f[i] = 0; } else { Material * const &mat = shapes[i]->material; if (mat->texture) colour.setVector(i, mat->texture->evaluate(P.getVector(i))); else colour.setVector(i, mat->colour); ambient_f[i] = mat->ambient; diffuse_f[i] = mat->diffuse; specular_f[i] = mat->specular; shininess_f[i] = mat->shininess; } // ambient acc = colour * ambient; Shape *shadow_shapes[4]; vector<Light*>::iterator light; for (light = lights.begin(); light != lights.end(); light++) { // direction vector to light VectorPacket L = VectorPacket((*light)->pos) - P; L.normalize(); const mfloat4 L_dot_N = dot(L, N); mfloat4 valid = mcmpgt(L_dot_N, mZero); // test if this light is occluded (sharp shadows) if ((*light)->cast_shadows) { const RayPacket shadow_rays = RayPacket(P, L); union { mfloat4 dists; float dists_f[4]; }; dists = mInf; top->packet_intersection(shapes, shadow_rays, dists_f, shadow_shapes); valid = mand(valid, mcmpeq(dists, mInf)); } const VectorPacket R = L - N * mmul(mTwo, L_dot_N); const mfloat4 R_dot_V = dot(R, V); // diffuse acc.selectiveAdd(valid, colour * VectorPacket((*light)->colour) * mmul(diffuse, L_dot_N)); // specular valid = mand(valid, mcmpgt(R_dot_V, mZero)); mfloat4 spec = mmul(mmul(specular, mset1((*light)->colour.r)), mfastpow(R_dot_V, shininess)); acc.selectiveAdd(valid, spec); } return acc;}#endifvoid Raytracer::lightScatter(const Ray &ray, const Shape *shape, int depth, const Vector &P, const Vector &normal, bool from_inside, Colour &col){ if (depth < max_depth) { Colour trans_col, refl_col; Float trans = shape->material->transmissivity; Float refl = shape->material->reflectivity; const Float cos_i = - dot(normal, ray.dir); // reflection if (refl > 0.01) { Vector newdir = ray.dir + 2.0f * cos_i * normal; Ray newray = Ray(P, newdir); refl_col = raytrace(newray, depth + 1, shape); } // refraction if (trans > 0.01) { Float n, n1, n2; if (from_inside) { n1 = shape->material->refract_index; n2 = 1.0; n = n1; } else { n1 = 1.0f; n2 = shape->material->refract_index; n = 1.0f / 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.0f / (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; Vector newdir = n * ray.dir + (n*cos_i - cos_t) * normal; Ray newray = Ray(P + 0.001f*newdir, newdir); trans_col = raytrace(newray, depth + 1, NULL); } } col = (1-refl-trans)*col + refl*refl_col + trans*trans_col; } // ambient occlusion if (ao_samples && !from_inside) { Float miss = 0; for (int i = 0; i < ao_samples; i++) { Vector 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(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; }}Colour Raytracer::raytrace(Ray &ray, int depth, const 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 { const Vector P = ray.o + ray.dir * nearest_distance; // point of intersection Vector 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; } // shading function Colour col = PhongShader(nearest_shape, P, normal, ray.dir); lightScatter(ray, nearest_shape, depth, P, normal, from_inside, col); return col; }}#ifndef NO_SIMDvoid Raytracer::raytracePacket(RayPacket &rays, Colour *results){ union { float nearest_distances[4]; mfloat4 m_nearest_distances; }; mfloat4 mask; Shape *nearest_shapes[4]; static const Shape *origin_shapes[4] = {NULL, NULL, NULL, NULL}; m_nearest_distances = mInf; top->packet_intersection(origin_shapes, rays, nearest_distances, nearest_shapes); mask = mcmpneq(m_nearest_distances, mInf); if (!mmovemask(mask)) { for (int i = 0; i < 4; i++) results[i] = bg_colour; return; } const VectorPacket P = rays.o + rays.dir * m_nearest_distances; // point of intersection VectorPacket normal; for (int i = 0; i < 4; i++) if (nearest_shapes[i] != NULL) { const Vector Pvecti = P.getVector(i); normal.setVector(i, nearest_shapes[i]->normal(Pvecti)); } // make shapes double sided mfloat4 from_inside = mcmpgt(dot(normal, rays.dir), mZero); normal.mx = mselect(from_inside, msub(mZero, normal.mx), normal.mx); normal.my = mselect(from_inside, msub(mZero, normal.my), normal.my); normal.mz = mselect(from_inside, msub(mZero, normal.mz), normal.mz); // shading function VectorPacket pres = PhongShader_packet(nearest_shapes, P, normal, rays.dir); //pres.mx = mselect(mask, pres.mx, mset1(bg_colour.r)); //pres.my = mselect(mask, pres.my, mset1(bg_colour.g)); //pres.mz = mselect(mask, pres.mz, mset1(bg_colour.b)); for (int i = 0; i < 4; i++) if (nearest_shapes[i] != NULL) { results[i] = pres.getVector(i); lightScatter(rays[i], nearest_shapes[i], 0, P.getVector(i), normal.getVector(i), (mmovemask(from_inside)>>i)&1, results[i]); } else results[i] = bg_colour;}#endif#ifdef MSVC__declspec(noreturn)#else__attribute__((noreturn))#endifvoid *Raytracer::raytrace_worker(void *d){ static const int my_queue_size = 256; Raytracer *rt = (Raytracer*)d; Sample my_queue[my_queue_size]; Colour my_colours[my_queue_size]; int my_count; Ray ray;#ifndef NO_SIMD RayPacket rays; const bool can_use_packets = (rt->use_packets && rt->sampler->packetableSamples());#endif for (;;) { pthread_mutex_lock(&rt->sample_queue_mutex); while (rt->sample_queue_count == 0) { if (rt->end_of_samples) { dbgmsg(4, "T thread [%d] exiting\n", pthread_self()); pthread_mutex_unlock(&rt->sample_queue_mutex); pthread_exit(NULL); } pthread_cond_wait(&rt->sample_queue_cond, &rt->sample_queue_mutex); } if (rt->sample_queue_count >= my_queue_size) my_count = my_queue_size; else my_count = rt->sample_queue_count; rt->sample_queue_count -= my_count; // copy samples to local queue if (rt->sample_queue_pos + my_count >= rt->sample_queue_size) { register int c = rt->sample_queue_size - rt->sample_queue_pos; memcpy(my_queue, rt->sample_queue + rt->sample_queue_pos, c*sizeof(Sample)); memcpy(my_queue + c, rt->sample_queue, (my_count - c)*sizeof(Sample)); rt->sample_queue_pos = my_count - c; } else { memcpy(my_queue, rt->sample_queue + rt->sample_queue_pos, my_count*sizeof(Sample)); rt->sample_queue_pos += my_count; } if (rt->sample_queue_count <= my_queue_size*2) pthread_cond_signal(&rt->worker_ready_cond); pthread_mutex_unlock(&rt->sample_queue_mutex); // do the work#ifndef NO_SIMD if (can_use_packets) { // packet ray tracing assert((my_count & 3) == 0); for (int i = 0; i < (my_count >> 2); i++) { rt->camera->makeRayPacket(my_queue + (i<<2), rays); rt->raytracePacket(rays, my_colours + (i<<2)); } } else#endif { // single ray tracing for (int i = 0; i < my_count; i++) { ray = rt->camera->makeRay(my_queue[i]); my_colours[i] = rt->raytrace(ray, 0, NULL); } } // save the results pthread_mutex_lock(&rt->sampler_mutex); for (int i = 0; i < my_count; i++) rt->sampler->saveSample(my_queue[i], my_colours[i]); pthread_mutex_unlock(&rt->sampler_mutex); }}void Raytracer::render(){ if (!sampler || !camera || !top) return; static const int my_count_max = 256; sample_queue_size = my_count_max*2*num_threads; sample_queue = new Sample [sample_queue_size]; sample_queue_count = 0; int my_count, my_pos; int sampnum = 0, sampdone; int phase = 1; bool more_samples; sampler->init(); // create workers dbgmsg(1, "* using %d threads\n", num_threads); pthread_t *threads = new pthread_t[num_threads]; dbgmsg(1, "* raytracing...\n"); while ( (sampnum = sampler->initSampleSet()) > 0 ) { my_pos = 0; sample_queue_pos = 0; sampdone = 0; end_of_samples = false; for (int t = 0; t < num_threads; t++) { int rc = pthread_create(&threads[t], NULL, raytrace_worker, (void*)this); if (rc) { dbgmsg(0, "\nE pthread_create unsuccessful, return code was %d\n", rc); exit(1); } } dbgmsg(2, "- phase %d: 0%% done", phase); pthread_mutex_lock(&sampler_mutex); for (;;) { my_count = 0; while ( more_samples = sampler->nextSample(&sample_queue[my_pos++]) ) { my_count++; if (my_pos >= sample_queue_size) my_pos = 0; if (my_count >= my_count_max) break; } if (!more_samples && !my_count) break; pthread_mutex_unlock(&sampler_mutex); pthread_mutex_lock(&sample_queue_mutex); sample_queue_count += my_count; // wait for workers if there is enough samples ready on queue while (sample_queue_count > (2*num_threads-1)*my_count_max) pthread_cond_wait(&worker_ready_cond, &sample_queue_mutex); pthread_cond_signal(&sample_queue_cond); pthread_mutex_unlock(&sample_queue_mutex); sampdone += my_count; dbgmsg(2, "\b\b\b\b\b\b\b\b%2d%% done", ((long)sampdone - sample_queue_count)*100/sampnum); pthread_mutex_lock(&sampler_mutex); if (!more_samples) break; } dbgmsg(2, "\b\b\b\b\b\b\b\b100%% done\n"); pthread_mutex_unlock(&sampler_mutex); // wait for workers dbgmsg(2, "- waiting for threads to finish\n"); pthread_mutex_lock(&sample_queue_mutex); end_of_samples = true; while (sample_queue_count) { pthread_cond_broadcast(&sample_queue_cond); pthread_mutex_unlock(&sample_queue_mutex); pthread_mutex_lock(&sample_queue_mutex); } pthread_cond_broadcast(&sample_queue_cond); pthread_mutex_unlock(&sample_queue_mutex); for (int t = 0; t < num_threads; t++) pthread_join(threads[t], NULL); phase ++; } delete[] threads; delete[] sample_queue;}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 = Inf;}